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Image  segmentation  and  registration  are  vital  to  image  processing,  image 
analysis,  and  computer  vision.  In  the  last  decade,  research  into  these  areas  has 
developed  at  a rapid  pace  by  using  various  mathematical  methods.  In  this  dis- 
sertation, variational  partial  differential  equations  are  proposed  to  solve  image 
segmentation  and  registration  problems. 

A new  variational  partial  differential  equation  based  level  set  method  for 
simultaneous  image  segmentation  and  non-rigid  registration  is  presented.  This 
technique  incorporates  both  prior  shape  and  intensity  information.  Since  global 
rigid  registration  has  limited  applicability  when  non-rigid  shapes  are  considered, 
a transformation  is  created  which  is  the  sum  of  a global  rigid  function  and  a local 
non-rigid  deformation.  This  model  is  tested  against  two-chamber  end-systolic 
endocardial  ultrasound  images  of  thirteen  human  patients.  The  experimental 
results  provide  preliminary  evidence  of  the  effectiveness  of  the  model  in  detecting 
the  boundaries  of  the  incompletely  resolved  objects  which  were  plagued  by  noise, 
dropout,  and  artifact.  Develop  algorithms  for  generating  the  mean  shape  has  a 
significant  role  in  image  segmentation,  since  the  mean  shape  can  be  used  as  a prior 
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shape  to  acquire  better  segmentation  results.  As  alternate  ways,  several  statistical 
algorithms  including  a shape  related  energy  function.  Self-Organizing  map  with 
Procrustes  methods,  and  Self-Organizing  map  with  a principal  component  analysis 
to  generate  mean  shape  and  clustering  are  also  presented. 

Region  based  image  segmentation  and  registration  models  which  contains 
three  interrelated  sections:  image  segmentation  using  a modified  Mumford-Shah 
technique,  region  based  segmentation  using  a prior  shape,  and  simultaneous 
segmentation  and  registration  using  Mumford-Shah  model  are  also  proposed. 

The  first  goal  is  to  develop  an  image  segmentation  technique  using  a modified 
Mumford-Shah  model.  A variational  region  intensity  based  image  segmentation 
model  is  presented.  The  boundary  of  the  given  image  is  extracted  using  a modified 
Mumford-Shah  segmentation  technique.  Even  though  the  two  phase  case  is 
performed  for  image  segmentation  during  the  numerical  experiments,  the  suggested 
model  can  be  applied  to  more  general  cases.  Another  goal  is  to  develop  a region 
based  image  segmentation  technique  using  prior  shape  information.  The  prior 
shape  information  is  extracted  by  using  a modified  Mumford-Shah  segmentation 
technique.  The  prior  shape  knowledge  supports  the  segmentation  process  for  a 
novel  image.  Finally,  a region  based  model  for  simultaneous  image  segmentation 
and  registration  is  also  presented.  The  purpose  of  the  model  is  to  segment  and 
register  given  images  simultaneously  utilizing  a modified  Mumford-Shah  technique 
and  region  intensity  values.  The  segmentation  is  obtained  by  minimizing  a modified 
Mumford-Shah  model.  A global  rigid  registration  is  assisted  by  the  segmentation 
information  and  region  intensity  values.  A segmentation  and  registration  process  is 
obtained  simultaneously  in  this  model.  In  addition,  the  model  can  also  be  applied 
to  the  case  of  non-rigid  registration.  The  numerical  experiments  of  the  proposed 
models  are  tested  against  synthetic  data  and  simulated  human  brain  MRI  images. 
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The  experimental  results  provide  preliminary  evidence  of  the  effectiveness  of  the 
proposed  models. 


CHAPTER  1 
INTRODUCTION 

Developing  algorithms  for  solving  image  segmentation  and  registration  prob- 
lems have  become  an  important  and  active  area  of  research.  These  methods  have 
been  applied  to  a wide  range  of  problems  in  computer  vision,  pattern  recogni- 
tion, object  recognition,  image  processing,  and  image  analysis.  There  have  been  a 
multitude  of  different  approaches  for  separately  solving  image  segmentation  and 
registration  problems.  One  goal  of  this  research  is  to  develop  a method  which 
solves  segmentation  and  registration  problems  simultaneously. 

Segmentation  techniques  have  been  developed  to  capture  the  boundary  of  an 
object  by  several  different  approaches;  edge-based  methods  using  active  contour 
models  developed  by  Caselles,  Kimmel,  Sapiro,  Kass,  Olver,  Sethian,  Vemuri,  Yezzi, 
et  al.  [7,  8,  33,  38,  47,  86],  which  rely  on  the  information  of  the  edges,  namely  high 
magnitude  of  the  image  gradient;  region-based  methods  proposed  by  Chan.  Vese, 
Cremers,  Tsai,  Mumford,  Shah,  Yezzi,  Zhu,  et  al.  [9,  10,  19,  22,  57,  75,  87,  93] 
that  depend  on  the  region  statistics  or  homogeneity;  or  a combination  of  the  two 
presented  by  Paragios  and  Deriche  [60,  61,  62]  using  Geodesic  Active  Region 
models.  However,  these  methods  may  not  be  able  to  accommodate  all  types 
of  imaging  difficulties  that  occur  in  images,  including  significant  signal  loss, 
noise,  and  non-uniformity  of  regional  intensities.  Therefore,  prior  information  is 
necessary  to  make  the  resulting  segmentation  more  accurate.  The  prior  shape  and 
intensity  information  has  been  incorporated  into  the  geodesic  active  contour  model 
presented  by  Chen,  Leventon,  Faugeras,  et  al.  [12,  13,  14,  15,  41,  42].  The  details 
are  available  in  Chen,  Huang,  et  al.  [13].  Recently,  the  prior  shape  and  intensity 
information  has  been  also  incorporated  into  deformable  models  proposed  by 


1 


2 


Metaxas,  Huang,  et  al.  [31,  32].  A class  of  deformable  models  called  MetaMorphs 
which  possess  both  shape  and  interior  texture  statistics  integrate  boundary  and 
region  information  coherently  in  a common  variational  framework  [32].  The 
MetaMorphs  models  [32]  are  extended  into  a unified  shape  and  intensity  space 
by  implicitly  representing  shapes  as  images  in  the  space  of  distance  transforms. 

A novel  stochastic  chord  based  matching  algorithms  align  global  transformation 
considering  both  shape  and  gray-level  intensity  information.  The  complementary 
local  registration  is  performed  by  deforming  a Free  Form  Deformation  control 
lattice  to  maximize  mutual  information  between  both  shape  and  intensity  images 

[31]- 

Image  registration  is  the  process  of  overlaying  two  or  more  images  of  the  same 
scene  or  similar  scene  taken  at  different  times,  and/or  by  different  sensors.  Area 
based  methods  or  feature  based  methods  have  been  developed  to  match  given 
images.  For  the  details  of  the  image  registration  techniques,  refer  to  [92]  presented 
by  Zitova  and  Flusser.  Multi-modal  registration  combined  with  mutual  information 
(MI)  proposed  by  Hata,  Maes,  Viola,  Wells,  et  al.  [27,  46,  80,  83]  generated  a new 
class  of  rigid/affine  registration  algorithms.  Nonlinear  registration  is  the  set  of 
techniques  that  allow  the  alignment  of  data  sets  that  are  mismatched  in  a nonlinear 
or  nonuniform  manner.  Dense  correspondences  between  pairs  and  groups  of  images 
were  generated  using  non-rigid  algorithms  presented  by  Bro-Nielsen,  Chefd'Hotel, 
Faugeras,  Gee,  Marsland.  Rueckert,  et  al.  [5,  11,  25,  51,  69].  Several  approaches 
have  been  made  in  establishing  the  application  of  mutual  information  for  non-rigid 
registration  proposed  by  Gaens,  Hermosillo,  Hellier,  Likar,  Meyer,  Rueckert,  et  al. 
[24,  28,  29,  44.  54,  68],  Further  improvements  of  these  methods  have  been  done  by 
either  the  inclusion  of  some  prior  information  on  the  joint  intensity  distribution 
(e.g.  Likar  and  Permus  [44])  or  the  use  of  higher  order  definitions  of  mutual 
information  which  incorporates  spatial  information  (e.g.  Rueckert,  Clarkson,  et 
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al.  [68]).  Instead  of  intensity  information,  a feature  based  model  combines  region 
matching  and  segmentation  with  free  boundary  conditions  in  Richard  and  Cohen 
[67].  However,  the  model  in  Richard  and  Cohen  [67]  can  be  applied  only  in  a single 
region  of  interest. 

However,  segmentation  or  registration  technique  alone,  dependence  of  seg- 
mentation on  registration  methods,  or  dependence  of  registration  on  segmenta- 
tion methods  does  not  solve  the  image  segmentation  and  registration  problems 
completely.  Recently,  joint  segmentation  and  registration  methods  have  been  de- 
veloped through  the  use  of  active  contours  proposed  by  Chen,  Yezzi,  Young,  et 
al.  [13,  15,  88,  89,  90],  An  explicit  combination  of  registration  with  segmentation 
has  been  developed  in  a variational  framework  through  active  contours  by  Yezzi, 
Zollei.  and  Kapur  [88,  89].  Moelich  and  Chan[56]  extend  the  algorithms  of  Yezzi, 
Zollei,  and  Kapur  [88]  to  joint  segmentation  and  of  an  object  in  two  images,  and 
the  prior  knowledge  is  implicitly  incorporated  through  the  simultaneous  evolution 
of  the  registration  maps  and  the  active  contour  in  the  variational  framework.  In 
Chen,  Tagare,  et  al.  [15],  the  shape  model  was  created  by  averaging  the  aligned 
training  samples.  Better  maps  were  gathered  by  matching  the  evolving  interface 
to  the  shape  prior.  The  geometric  active  contour  was  minimized  and  the  shape 
disparity  energy  function  between  the  active  contour  and  prior  shape  was  measured 
by  a distance  function.  This  algorithm  was  extended  in  Chen,  Guo,  et  al.  [14]  by 
incorporating  the  information  on  the  location  of  a few  key  points  gathered  from 
standardized  results.  Instead  of  using  a few  key  points  by  Chen.  Guo,  et  al.  [14], 
an  intensity  profiling  term  was  added  to  match  the  intensity  of  the  interface  to 
the  prior  intensity  cross  prior  shape  and  the  active  contour  in  Chen,  Huang,  et  al. 
[13].  The  mutual  information  intensity  profiling  and  the  determination  of  the  pa- 
rameter which  balances  the  influence  from  image  information  and  shape  prior  were 
developed  in  Chen,  Huang,  et  al.  [12]  similar  to  Chen,  Huang,  Wilson,  et  al.  [13]. 
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The  morphing  active  contours  algorithm  is  combined  with  the  joint  segmentation 
and  registration  to  deal  with  the  complexity  of  segmenting  real  images  with  an 
application  to  CT  scans  for  radiotherapy  treatment  planning  by  Young  and  Levy 

[90]- 

The  purpose  of  Chapter  2 is  to  develop  an  algorithm  to  solve  simultaneous 
image  segmentation  and  registration  problem.  An  edge-based  active  contour  model 
using  prior  shape  and  intensity  profile  is  suggested  for  solving  simultaneous  seg- 
mentation and  registration  problem  in  Chapter  2.  In  Chapter  2,  a new  variational 
PDE  based  level  set  method  for  simultaneous  image  segmentation  and  registration 
is  presented  with  an  application  to  two-chamber  human-heart  end-systolic  endo- 
cardial ultrasound  images.  The  purpose  of  the  model  is  to  segment  the  endocardial 
borders  of  the  given  images  using  prior  shape  and  intensity  profile.  Prior  shape  and 
prior  intensity  of  thirteen  human-heart  two-chamber  epicardial  ultrasound  images 
were  generated  by  the  method  in  Chen,  Feng,  et  al.  [12],  The  segmentation  is 
obtained  by  finding  a non-rigid  registration  to  the  prior  shape.  The  model  performs 
a simultaneous  segmentation  and  registration  in  a way  similar  to  Chen,  Yezzi,  et 
al.  [13,  15,  88]  through  active  contours.  But  the  model  differs  from  Chen,  Yezzi, 
et  al.  [13,  15,  88]  by  using  a shape  prior  as  an  initial  contour.  Since  the  shape 
prior  is  used  as  an  initial  contour,  the  level  set  form  of  the  model  is  fixed  which 
does  not  require  the  re-initialization  process  during  the  numerical  calculation.  The 
model  in  this  chapter  is  similar  to  Chen,  Leventon,  et  al.  [12,  13,  41]  , because  it 
uses  both  prior  shape  and  prior  intensity  information  to  produce  the  segmentation. 
The  segmentation  in  this  model  is  obtained  by  a non-rigid  deformation  from  the 
prior  shape.  This  model  follows  the  idea  from  Huang.  Metaxas,  Paragios,  Yezzi, 
et  al.  [31,  32,  63.  64,  71]  for  matching  nonequivalent  shapes  by  the  combination 
of  a global  rigid  transformation  and  a local  non-rigid  deformation.  However,  the 
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model  differs  from  Paragios.  Rousson,  and  Ramesh  [63,  64]  by  having  an  inten- 
sity matching  term  in  the  energy  function  which  provides  the  difference  of  the 
intensity  variation  across  the  prior  shape  and  the  shape  of  interest.  Minimizing 
the  difference  of  the  gradients  of  the  intensity  between  a smoothed  prior  image 
and  a smoothed  given  image  is  more  accurate  than  minimizing  the  gradient  term 
in  the  given  image,  due  to  the  loss  of  the  information  of  the  gradient  near  the 
edge  of  the  desired  object.  Generating  the  mean  shape  has  an  important  appli- 
cation to  image  segmentation,  since  it  can  be  used  as  a prior  shape  to  overcome 
all  types  of  imaging  difficulties.  Chen,  leventon,  Paragios,  Yezzi,  Duncan,  et  al. 

[12,  13,  14,  15,  41,  42,  63,  64,  71,  85]  showed  that  prior  mean  shape  can  play  an  im- 
portant role  to  acquire  better  segmentation  results.  Several  statistical  methods  for 
generating  the  mean  shape  and  clustering  are  also  presented  as  alternate  techniques 
in  Chapter  2.  To  create  the  mean  shape,  the  shape  related  energy  function  is  mini- 
mized with  an  orthogonal  Procrustes  measurement.  The  other  model  is  introduced 
to  generate  the  mean  shape  and  clustering  simultaneously  using  a Self-Organizing 
map  combined  with  Procrustes  methods.  To  measure  the  similarity  between  two 
contours,  an  area  difference  distance  measurement  proposed  by  Chen,  Wilson, 
and  Feng  [16]  and  orthogonal  Procrustes  measurement  [73]  are  used.  In  addition, 
a Self-Organizing  map  combined  with  a Principal  Component  Analysis  model  is 
also  introduced  to  generate  a better  mean  shape.  These  algorithms  are  applied 
to  human  heart  epicardial  cardiac  borders  trace  by  an  expert  on  two-chamber 
end-systolic  echo-cardiographic  images.  For  verifying  these  results,  mean  distance, 
standard  deviation,  and  Tanimoto  distance  measurement  [40]  are  used  to  compare 
the  results. 

Chapter  2 is  organized  as  follows:  active  contours  are  reviewed  briefly  in 
section  2.2.  Generating  prior  shape  and  intensity  information  are  explained  in 
section  2.3.  In  section  2.4,  the  model  is  proposed  for  a simultaneous  segmentation 
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and  non-rigid  registration  and  the  level  set  form  of  the  proposed  model  is  described 
in  section  2.5.  In  section  2.6.  numerical  methods  pertaining  to  the  model  are 
explained.  Experimental  results  of  the  model  which  were  applied  to  human- 
heart  two-chamber  end-systolic  ultrasound  images  are  also  shown  in  section  2.7. 
Several  statistical  methods  with  an  application  to  human  heart  cardiac  borders  for 
generating  the  mean  shape  and  clustering  are  presented  in  section  2.8. 

The  goal  of  Chapter  3 is  to  formulate  various  region  based  algorithms  for 
solving  image  segmentation  and  registration  problems.  Several  region-based 
algorithms  for  image  segmentation  and  registration  are  proposed  in  Chapter  3, 
which  contains  three  interrelated  sections:  image  segmentation  using  a modified 
Mumford-Shah  technique,  region  based  image  segmentation  with  prior  shape,  and  a 
modified  Mumford-Shah  model  based  simultaneous  segmentation  and  registration. 

Mumford  and  Shah  introduced  the  most  celebrated  region  based  image 
segmentation  model  in  1989  [57].  The  measurement  of  an  edge  length  term  in  the 
Mumford-Shah  model  is  approximated  by  a quadratic  integral  of  an  edge  signature 
function  proposed  by  Ambrosio  and  Tortorelli  in  1990  [1].  Chan  and  Vese  proposed 
a piecewise  constant  Mumford-Shah  model  in  [9,  10]  Developments  of  variational 
level  set  implementation  techniques  based  Mumford-Shah  model  are  followed  by 
Esedoglu,  Gibou,  Likar,  et  al.  [22,  26,  43].  The  segmentation  is  represented  by 
characteristic  functions  using  phase  fields  in  Esedoglu,  Likar,  et  al.  [22,  43].  The 
details  of  phase  field  theory  can  be  found  in  Baldo,  Modica,  Shen,  Zhou,  et  al. 

[3,  55,  70,  82]. 

Recently,  a soft  Mumford-Shah  segmentation  is  introduced  in  Shen  [70].  The 
model  in  this  section  is  motivated  by  Esedoglu,  Shen,  et  al.  [22,  70,  77]. 

In  section  3.1,  a region  based  variational  partial  differential  equation  model 
for  an  image  segmentation  is  presented  with  an  application  to  synthetic  image 
and  simulated  human  brain  MRI  images.  The  purpose  of  the  model  is  to  segment 
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borders  of  the  given  images.  The  segmentation  is  obtained  by  using  a modified 
Mumford-Shah  segmentation  technique.  In  section  3.1.2,  a modified  Mumford-Shah 
model  is  proposed.  The  Euler-Lagrange  equations  of  the  suggested  model  are 
presented  in  this  section  as  well.  In  section  3.1.3,  numerical  methods  pertaining 
the  proposed  model  are  explained.  The  model  is  applied  to  synthetic  data  and 
simulated  human  brain  MRI  images  are  shown  in  section  3.1.4. 

In  section  3.2,  a new  variational  image  segmentation  model  based  on  the  region 
intensity  information  utilizing  shape  knowledge  is  suggested.  The  goal  of  the  model 
is  to  segment  a given  novel  image  using  prior  shape  information.  The  prior  shape 
knowledge  is  obtained  by  using  modified  Mumford-Shah  segmentation  technique  in 
section  3.1.  The  prior  shape  knowledge  is  assisting  the  image  segmentation  process 
to  a given  novel  image.  In  section  3.2.2.  the  model  for  generating  prior  shape 
information  using  a modified  Mumford-Shah  segmentation  technique  is  described. 

In  section  3.2.3,  the  proposed  model  is  introduced.  The  Euler-Lagrange  equations 
of  the  suggested  model  are  presented  in  this  section  as  well.  The  model  is  applied 
to  synthetic  data  and  simulated  human  brain  MRI  images  are  shown  in  section 
3.2.4. 

In  section  3.3,  a modified  Mumford-Shah  model  based  simultaneous  segmen- 
tation and  registration  technique  is  presented.  The  purpose  of  the  model  is  to 
segment  and  register  given  images  simultaneously  utilizing  modified  Mumford-Shah 
technique  and  region  intensity  values.  The  model  performs  a simultaneous  segmen- 
tation and  registration  in  a similar  way  from  Yezzi,  Zollei,  and  Kapur  [88,  89].  But 
the  model  differs  from  Yezzi,  Zollei,  and  Kapur  [88,  89]  by  using  a region  intensity 
and  modified  Mumford-Shah  segmentation  technique.  In  section  3.3.2,  the  model  is 
proposed  for  a simultaneous  segmentation  and  registration  and  the  Euler-Lagrange 
equations  of  the  proposed  model  are  described.  In  section  3.3.3,  numerical  methods 
pertaining  to  the  model  are  explained.  The  model  is  applied  to  synthetic  data  and 
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simulated  human  brain  MRI  images  are  shown  in  section  3.3.4.  Finally  conclusions 
and  future  work  are  stated  in  Chapter  4. 


CHAPTER  2 

AN  EDGE-BASED  ACTIVE  CONTOUR  MODEL  WITH  PRIOR  SHAPE  AND 

INTENSITY  PROFILE  FOR  SIMULTANEOUS  SEGMENTATION  AND 
NON-RIGID  REGISTRATION 

2.1  Introduction 

The  goal  of  image  segmentation  is  to  capture  the  boundary  of  an  object. 

Using  the  information  of  the  edges,  namely  high  magnitude  of  the  image  gradient, 
edge-based  methods  have  been  developed  using  active  contour  models  proposed 
by  Caselles,  Kimmel,  Sapiro,  Kass,  Olver,  Sethian,  Vemuri,  Yezzi,  et  al.  [7,  8, 

33,  38,  47,  48,  50,  53,  86].  To  accommodate  all  types  of  imaging  difficulties  that 
occur  in  cardiac  ultrasound  images  including  significant  signal  loss,  noise,  and 
non-uniformity  of  regional  intensities,  the  prior  shape  and  intensity  information 
has  been  incorporated  into  the  geodesic  active  contour  model  presented  by  Chen, 
Leventon,  Faugeras,  et  al.  [12,  13,  14,  15,  41,  42],  Recently,  the  prior  shape  and 
intensity  information  has  been  also  incorporated  into  deformable  models  proposed 
by  Metaxas.  Huang,  et  al.  [31,  32], 

Image  analysis  has  been  developed  rapidly  with  applications  to  medical 
imaging,  target  tracking,  and  computer  vision.  Especially  in  medical  imaging  and 
target  tracking,  image  registration  play  a significant  role.  Area  based  methods  or 
feature  based  methods  have  been  developed  for  image  registration.  Hata,  Maes, 
Viola,  Wells,  et  al.  [27,  46,  80,  83]  proposed  multi-modal  registration  combined 
with  mutual  information  (MI).  This  brought  a new  class  of  rigid/affine  registration 
algorithms.  Nonlinear  registration  is  the  set  of  techniques  that  allow  the  alignment 
of  data  sets  that  are  mismatched  in  a nonlinear  or  nonuniform  manner.  Several 
approaches  have  been  made  in  establishing  the  application  of  mutual  information 
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for  non-rigid  registration  presented  by  Gaens,  Hermosillo.  Hellier,  Likar,  Meyer, 
Rueckert,  et  al.  [24,  28,  29,  44,  54,  68].  For  the  details  of  the  image  registration 
techniques,  please  refer  to  [92]  presented  by  Zitova  and  Flusser. 

Recently,  joint  segmentation  and  registration  methods  have  been  developed 
to  get  better  image  segmentation  and  registration  results  through  active  contours 
proposed  by  Chen,  Yezzi,  Young,  et  al.  [13,  15,  88,  89,  90]. 

In  Chapter  2,  a new  variational  partial  differential  equation  (PDE)  based  level 
set  method  for  simultaneous  image  segmentation  and  registration  is  presented  with 
an  application  to  two-chamber  human-heart  end-systolic  ultrasound  images.  The 
purpose  of  the  model  is  to  segment  the  endocardial  borders  of  the  given  images 
using  prior  shape  and  intensity  information.  Prior  shape  and  prior  intensity  of 
thirteen  human-heart  two-chamber  end-systolic  ultrasound  images  were  generated 
by  the  method  in  Chen,  Feng,  et  al.  [12].  The  segmentation  is  obtained  by  finding 
a non-rigid  registration  to  the  prior  shape.  The  model  performs  a simultaneous 
segmentation  and  registration  using  an  active  contour  approach  similar  to  the  one 
used  by  Chen,  Yezzi,  et  al.  [13,  15,  88].  But  the  model  differs  from  Chen,  Yezzi  et 
al.  [13,  15,  88]  by  using  a shape  prior  as  an  initial  contour.  Since  the  shape  prior  is 
used  as  an  initial  contour,  the  level  set  form  of  the  model  is  fixed  which  does  not 
require  the  re-initialization  process  during  the  numerical  calculation.  The  model  in 
this  Chapter  is  similar  to  the  one  used  by  Chen,  Leventon,  et  al.  [12,  13,  41],  since 
it  uses  both  prior  shape  and  prior  intensity  information.  However  the  segmentation 
in  our  model  is  obtained  by  finding  the  non-rigid  deformation  to  the  prior  shape. 
The  model  in  this  Chapter  follows  the  idea  from  Paragios,  Metaxas,  Yezzi,  et.  al. 
[31,  32,  63,  64,  71]  for  matching  nonequivalent  shapes  by  the  combination  of  a 
global  rigid  transformation  and  a local  non-rigid  deformation.  However,  the  model 
differs  from  Paragios  [63,  64]  by  having  an  intensity  matching  term  in  the  energy 
function  which  provides  the  difference  of  the  intensity  variation  across  the  prior 
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shape  and  the  shape  of  interest.  Minimizing  the  difference  of  the  gradients  of  the 


intensity  between  a smoothed  prior  image  and  a smoothed  given  image  is  more 
accurate  than  minimizing  the  gradient  term  in  the  given  image,  due  to  the  loss  of 
the  information  of  the  gradient  near  the  edge  of  the  desired  object. 

Chapter  2 is  organized  as  follows:  active  contour  is  reviewed  briefly  in  section 
2.2.  Generate  prior  shape  and  intensity  information  is  explained  in  section  2.3.  In 
section  2.4,  the  model  is  proposed  for  a simultaneous  segmentation  and  non-rigid 
registration.  The  level  set  form  of  the  proposed  model  and  Euler-Lagrange  equa- 
tions of  the  level  set  form  are  described  in  section  2.5.  In  section  2.6,  numerical 
methods  pertaining  to  the  model  are  explained.  Experimental  results  of  the  model 
which  were  applied  to  human-heart  two-chamber  end-diastolic  ultrasound  images 
are  also  shown  in  section  2.7.  Alternate  ways  to  create  prior  mean  shape  and 
clustering  using  various  statistical  methods  are  shown  in  section  2.8. 


Active  contours  are  used  to  control  the  smoothness  of  the  curve  and  to  attract 
the  contour  towards  the  boundary  of  the  object,  especially  they  are  used  for 
image  segmentation  which  is  finding  object  boundaries  in  images.  In  this  section, 
edge-based  active  contour  models  are  briefly  reviewed: 


be  an  image  defined  on  Q.  Let  C(p)  = (x(p),y(p))(p  £ [0, 1])  be  a differentiable 
parameterized  curve  in  12.  Kass,  et  al.  [33]  proposed  the  classical  snake  model 
which  is  formulated  by  minimizing  the  energy  function: 


The  first  term  in  Equation  (2-1)  is  called  the  internal  energy  and  is  used  to 


2.2  Background  Review;  Active  Contours 


Let  £l  in  R 2 be  a closed  and  bounded  region  of  the  plane  and  I : —>  R 


(2-1) 


control  the  smoothness  of  the  curve.  The  second  term  is  the  external  energy  and  is 
used  to  attract  the  contour  towards  the  boundary  of  the  object. 
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The  following  geometric  deformable  contour  model  based  on  the  level  set 
interpretation  of  the  Euclidean  curve  was  proposed  by  Caselles,  Malladi,  et  al. 
[7,  8,  47]: 


ut  = 5(|V/|)|Vu|(dfu(— ) + u),  (2-2) 

where  y(|V/|)  can  be  chosen  as 

= 1 + /3\S7(Ga  * I)\2  ’ 

|a;|2 

where  (3  > 0 and  u > 0 are  parameters,  Ga{x)  = ^ e~  , so  that  <?(| V/|)  is  strictly 
positive  in  homogeneous  regions  and  near  zero  on  the  edges  proposed  by  Caselles, 
Cootes,  Yezzi,  et  al.  [8,  17,  86].  In  this  formation,  the  evolving  contour  is  the  zero 
level  set  ( u(x,y,t ) = 0)  of  the  solution  of  Equation  (2-2).  The  level  set  formulation 
is  an  improvement  on  the  original  because  it  can  handle  multiple  contours  and 
topological  changes. 

A further  modification  called  geometric  active  contour  model  was  proposed  by 
Caselles,  Kichenassamy,  et  al.  [8,  38].  The  main  modification  was  to  propose  an 
image  gradient  dependent  Riemannian  arc  length  function  for  the  contour.  The  arc 
length  was  taken  to  be  the  ’’energy”  E(C ) of  the  contour  and  was  designed  in  such 
a way  that  minimizing  it  leads  to  active  contour  becoming  stationary  at  pixels  of 
high  gradient: 

E(C)  = f 5(|V/|)(C(p))|C"(p)|dp,  (2-3) 

Jo 

where  g(|V/|)  is  the  function  defined  in  Equation  (2-3).  Because  of  the  change  of 
arc-length  measure  term  \C'(p)\dp,  the  model  does  not  depend  on  the  choice  of  the 
parameters  any  more. 
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2.3  Generating  Prior  Shape  and  Intensity  information 

The  purpose  of  our  suggested  model  is  to  capture  the  boundary  of  the  given 
images  by  using  the  prior  shape  and  prior  intensity  information.  It  has  been  shown 
that  prior  mean  shape  can  play  an  important  role  to  acquire  better  segmentation 
results  presented  by  Chen,  leventon,  Paragios,  Yezzi,  Duncan,  et  al.  [12,  13,  14, 

15,  41,  42,  63,  64,  71,  85].  Prior  shape  and  prior  intensity  in  our  model  were 
generated  by  the  method  in  Chen,  Feng,  et  al.  [12]  with  an  application  to  thirteen 
human-heart  two-chamber  epicardium  ultrasound  images.  The  followings  are  brief 
descriptions  of  the  method  to  generate  the  prior  mean  shape  and  clustering  from 
Chen,  Feng,  et  al.  [12]  : 

2.3.1  Area  Difference  Distance  Measurement 

There  are  several  different  ways  to  define  the  distance  between  two  contours. 
The  area  difference  distance  measurement  is  used  from  Chen,  et  al.  [16]  to  get 
the  distance  between  two  contours.  The  following  is  the  brief  review  of  the  area 
distance  and  the  model  developed  in  Chen,  et  al.  [16]: 

The  notion  of  shape  in  our  model  is  assumed  independent  of  translation, 
rotation,  and  scaling. 

Definition  1.  Let  A and  B be  two  contours.  Then  area(A,  B ) = the  area  of 
(Int(A)  U Int(B)  — Int(A)  fl  Int(B)),  where  the  ’’area”  of  digitized  curves  is  defined 
and  computed  as  in  Chen  et  al.  [16]. 

Definition  2.  The  area  distance  between  two  contours  A and  B is  defined  as 
the  minimum  area  between  two  contours  after  rotation,  translation,  and  scaling. 

AD(A,B):—  min  (area(nRA  + T,  B))  (2-4) 

a —b 


In  this  subsection,  we  will  denote  fiR  = R = 


b a 
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Theorem  1.  Let  A = (a^),  5 = (&„),  and  ct  — (bi+ 1,2  - bia)2  + (&i+i,i  - &i,i)2, 
where  i = 1 n and  j = 1,2.  Then 


mm(area(RA  + T,  B )) 

R,T 


is  achieved  at  R — 


a —b 
b a 


, T = 


t i 
t2 


, where 


a 

51 

b 

52 

53 

^2 

54 

Cl  0 C 2 C3 
0 Cl  — C3  C2 
C2  — C3  C4  0 
C3  C2  0 C4 


Cl  = ^2ci(a2l  + a22),  C2  = ^CjOii,  C3  = ^Ciai2,  C4  = ^Cj, 

i=l  i=l  t=l 

n n 

51  = Cj(6jiaii  + b^att),  B 2 = — bua^), 

1=1 

n n 

B3  = Cibn , 54  = Cifei2. 


(2-5) 


i=i 


i=i 


i=l  i=l 

Here,  area(RA  + T,B)  was  approximated  by 


'y  ^ Ct((aa,x  — bai2  + t\  — 5ji)2  + (6aa  + aflt2  + t2  — bj2)2). 


1=1 


2.3.2  Generating  Prior  Mean  Shape  and  Clustering 

The  method  to  generate  prior  mean  shape  and  clustering  proposed  by  Chen, 
et  al.  [16]  is  briefly  reviewed.  To  get  the  prior  mean  shape,  first  choose  the  first 
one  of  given  n contours  as  an  initial  mean  shape.  Then  align  the  rest  of  contours 
to  an  initial  contour  to  obtain  Rx,. . . , Rn,  and  Tx, ...  ,Tn  using  the  area  difference 
distance  measurement  from  section  2.3.1.  Finally,  the  prior  mean  shape  C*  is 
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created  by 


n 


C*  = (Cl  + Y/RjCj  + TJ)/n. 


(2~6) 


If  given  contours  have  large  shape  variations,  it  is  necessary  to  have  a clustering 
process  to  obtain  better  prior  mean  shape.  If  the  n training  contours  Cj,  j — 

1, . . . ,n  are  to  be  grouped  into  k clusters,  take  k arbitrary  contours  as  the  initial 
contours.  For  the  first  cluster,  fix  an  arbitrary  contour,  say  C\.  Then  align  the  rest 
of  the  contours  to  C\  using  Equation  (2-4)  to  obtain  RjCj  + Tj,  correspondingly 
AD(Ci,Cj),j  = l,...,n  — l.  Define 


Then  the  first  cluster  consists  of  all  the  curves  Ck,  such  that  AD(C\,  Ck)  < cr/3. 
The  remaining  contours  are  grouped  into  k — 1 clusters  by  repeating  the  above 
process.  A prior  intensity  information  image  is  created  by  using  mutual  information 
of  image  geometry.  The  details  can  be  found  in  Chen,  Huang,  et  al.  [12].  These 
algorithms  are  tested  against  eighty-five  two-chamber  end-systolic  endocardial 
human-heart  cardiac  borders.  These  contours  are  grouped  into  three  clusters.  The 
algorithms  to  generate  prior  mean  shape  and  intensity  information  is  tested  against 
one  of  the  clusters  which  has  thirteen  contours.  Prior  mean  shape  and  intensity 
information  of  these  thirteen  contours  and  their  corresponding  images  are  created 
using  the  method  from  Chen,  Feng,  et  al.  [12].  Prior  mean  shape  C*  and  intensity 
information  I*  of  these  thirteen  contours  and  their  corresponding  images  which 
have  0.62 mm  resolution  in  each  pixel  are  created.  The  prior  intensity  information 
and  a prior  shape  as  a solid  line  are  in  Figure  2-1. 


The  purpose  of  our  model  is  to  capture  the  boundary  of  the  given  image  by 
using  the  prior  shape  and  prior  intensity  information.  Our  prior  shape  and  prior 
intensity  of  thirteen  human-heart  two-chamber  epicardial  ultrasound  images  were 


n 


(2-7) 


2.4  Description  of  the  Proposed  Model 
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generated  by  the  method  from  Chen  , Huang  et  al.  [12].  The  model  is  proposed  as 
follows: 

Find  a,  u,  v,  y,  R,  T by  minimizing  the  energy  function 

E(a,  u,  v,  y,  R,  T)  = \\  fN(c.}  | V(Ga*  I*)(x,  y)  - aV(Ga  * I)(u,  v)\2dx 

+A2  /jv(c*)(M2  + \v\2)dx  + A3  /JV(c.)(|Vti|2  + | Vv| 2)dx  (2-8) 

+A4  fo  5(|V/|)(ft, v)(C*(p))\(C*)'(p)\dp, 

5(|V/|)  can  be  chosen  as  5(|V/|)  = » P > 0 is  a parameter, 

u(x,  y) 


Ga(x)  = and 


u{x,y) 
v (x,y) 


yR 


y 


+ T + 


v(x,y) 


An  average  image 


I*,  an  original  image  I,  and  the  average  curve  C*  axe  given  and  Ai,A2,  A3,  and  A4 
are  parameters.  N(C*)  is  the  neighborhood  of  the  prior  shape  C*  and  a,  y,  R,  and 
T are  constants. 

The  first  term  minimizes  the  difference  of  the  intensity  variation  across  the 
prior  shape  and  the  shape  of  interest.  Minimizing  the  difference  of  gradients  of  the 
intensity  between  smoothed  prior  image  and  the  smoothed  given  image  is  more 
accurate  than  minimizing  the  gradient  term  in  the  given  image,  due  to  the  loss  of 
the  information  of  the  gradient  near  the  edge  of  the  desired  object.  Since  images 
may  have  different  levels  of  intensity,  the  model  is  also  minimized  over  the  scale 
factor  variable  a.  The  second  term  is  minimizing  the  magnitude  of  the  non-rigid 
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deformation  term 


. Smoothing  the  non-rigid  term 


is  in  the  third  term. 


Final  term  is  the  energy  functional  of  active  contours  proposed  by  Caselles,  Yezzi, 
et  al.  [8,  38],  that  measures  the  amount  of  high  gradient  under  the  trace  of  the 
curve  and  increases  the  attraction  of  the  evolving  curve  toward  the  boundary  of 
the  object.  Since  the  prior  shape  C*  is  used  as  an  initial  contour  and  the  domain 
is  taken  as  the  neighborhood  of  the  prior  shape,  it  helps  the  numerical  solutions 
converge  faster.  In  addition,  when  the  model  is  formulated  into  the  level  set  form  in 
Equation  (2-9),  the  level  set  of  the  active  contour  in  this  model  is  fixed.  Therefore, 
the  re-initializing  process  in  the  calculation  is  not  needed  which  decreases  the 
numerical  calculation  time.  In  this  model,  a simultaneous  segmentation  and 
non-rigid  registration  is  performed.  The  segmentation  is  obtained  by  finding  a 
non-rigid  registration  to  the  prior  shape.  A registration  consists  of  a global  rigid 
transformation  and  a local  non-rigid  deformation. 

2.5  Level  Set  Formulation  of  the  Proposed  Model 
2.5.1  Level  Set  Formulation 

The  level  set  methods  proposed  by  Osher,  Sethian,  et  al.  [59]  are  extensively 
used  in  the  problems  of  curve  evolution,  because  they  allow  the  curve  to  develop 
cusps,  corners,  and  topological  changes.  Now  we  give  a variational  level  set 
approach  proposed  by  Chan,  Vese,  Osher,  et  al.  [9,  91]  : 

Let  the  contour  C be  the  zero  level  set  of  a Lipschitz  function  w such  that 
{j:|u;(x)  > 0}  is  the  set  inside  C.  Let  H(z ) be  the  Heaviside  function,  that  is 
H(z ) = 1 if  z > 0,  and  H(z ) = 0 if  z < 0,  and  <5  = H'(z)( in  the  sense  of 
distribution)  be  the  Dirac  measure.  Then,  the  length  of  the  zero  level  set  of  u 
in  the  conformal  metric  ds  = g\C'(p)\dp  can  be  measured  by  fn  H(w)\)  = 

where  fl  is  the  image  domain.  Therefore,  the  level  set  formulation 
of  our  variational  approach  is  to  find  a,  u,  v,  //,  R,  T by  minimizing  the  energy 
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function 

E(a,u,v,n,R,T ) = Xi  Jfffc.)  |V(Gff  * I*)(x,y)  - aV(GCT  * I)(u,v)\2dx 


+a2  X 


AT(C-) 


it 


+ M2)dz  + A3  Xv(c-)(r Vu|2  + I Vv| 2)dx 


(2-9) 


+A4  Xv(C.,)  <5(c/c-)^(|  V/|)(?i,  -C)|  Vrfo*  |^, 

where  dc * is  the  distance  function  of  the  average  curve  G*,  g(|V/|)  can  be  cho- 

jj.  |2 

sen  as  g(|V/|)  = 1+ia[v(G  «/)p>  @ > 0 is  a parameter,  Ga(x ) = Xe-  and 

u(x.  11) 

. An  average  image  /*,  an  original  image 

I,  and  the  average  curve  C*  are  given  and  A x , A2 , A3,  and  A4  are  parameters. 

2.5.2  Euler-Lagrange  Equations: 

The  evolution  equations  associated  with  the  Euler-Lagrange  equations  of  (2-9) 
are 


u(x,y) 

--  y,R 

x 

+ T + 

u(x,  y) 

v(x,y) 

y 

v(x,y) 

* r)(x,  y ) - aV(Ga  * I)(u,  t)))(V(G<r  * I)(u,  v))dx 


(2-10) 


du 

dt 


= A!(V(Gct  * r)(x,y)  - aV(Ga  * I)(u,v)))a( 


d(V(Ga*I)(u,v)) 

du 


A 2u  + A3A u - A4J(dc. ) ( V ^ ) | Vrfc- 1 , in  N(C *) 


5n 

<9n 


du 

= 0,  on  dN(C*) 


(2-11) 


du 

dt 


= A1(V(Gff  * /*)(*,  y)  - aV(Ga  * 7)(fi,  f))))a(--(-V-(-^  * /)(?2’  5)) )- 

ov 


X2v  + A3Ao  - A4^(dc-)( 6,9(1  VaP(U’  U))[Vdc»|,  in  N{C *) 


do 

dn 


55 


= 0,  on  dN(C*) 


(2-12) 
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^ = A ! / (V(G,  * r)(x,  y)  - aV(G,  * /)(«, 

ot  Jn(c *) 


WWM)  + 3(V(G,»j)(i j,e» 


3u 


dv 


v)) 


dx- 


xj 

Jn 


IN{C ♦) 


3// 


|Vdc>|dx 


^ = Ai  f (V(G,  * /*)(*,  y)  - aV(Gff  * /)(«,  C)) 
at  Jn(c •) 


3(V(Gff*/)(u,u))  3(V(Gff  */)(«,(;))  dfl 
3u  3u  jdd 


a: 

y 


ydx— 


xj 

Jn(C *) 


3/i 


dd 


/x|Vdc*|dx 


3T 

3t 


Ax  [ ( V(G,  * /*)(*,  y)  - aV(GCT  * /)(«,  0)) 

Jn(c*) 


,3(V(G,  */)(«, 0))  , 3(V(Gff  */)(«, 5))^, 

a( — 1 — J dx— 


du 


dv 


' N(C') 


3T 


(2-13) 


(2-14) 


(2-15) 


where  R is  the  rotation  matrix  in  terms  of  the  angle  9. 

2.6  Numerical  Methods 

In  this  section,  we  explain  how  to  solve  the  equations  (2-10)-(2-15)  numeri- 
cally. Equation  (2-9)  was  solved  by  finding  a steady  state  solution  of  (2-10)-(2-15). 
Following  the  approach  presented  by  Chan,  Vese,  et  al.  [9],  we  replace  H and  6 in 
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(2-10)-(2-15)  by  the  regularized  versions  of  them,  denoted 


He(z)  = 


and, 


5£(z)  = H'e{z) 


l,  if  z>  e 
0 , if  z < — e 

K1  + f + if  \z\  < e, 

0,  if  \z\  > e 

i[l  + CO«(t)]i  */  1*1  - £’ 


To  discretize  the  equations  (2-10)-(2-15),  we  use  a finite  difference  scheme 
and  we  applied  the  gradient  descent  method.  Let  iiy  denote  the  value  of  u at  the 
pixel  Xi  = ih  and  yj  = jk,  where  h and  k is  the  pixel  size.  Denote  u(xi,yj,tn)  by 
Uy.  The  time  derivative  ut  at  {i,j,tn)  is  approximated  by  the  forward  difference 
scheme: 

u"+i  _ un 

ut(i,j,tn)  = 

where  At  is  the  time  step.  Moreover,  (2- 10)-(2- 15)  were  implemented  by  the 
following  iteration  scheme  : 

u%,,u  — u”  u", . , n — Uii 

( \n  _ U+QJ  O \n  _ »U  + 1)  O 

\uxJii  ~ . • \uyJij  ~ , ) 


(«**)«  = 


U 


2u"  + it" 


h2 


) iuyy)ij  ~ 


<0+0  ~ 2u?-  + ufc.p 
k2 


((G„  * /)x)«  = 


(G„  * /)& 


(*+i)i 


(G<r  * /)"• 


lJ  h 

(G.  * /)»w+1)  - (Go-  * /)" 


((G„  * /),)£  = 
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(<?(|v/|U 


»(|V/|)?y+1)-»(|V/|)5 


y 


2.7  Experimental  Results 

It  is  very  difficult  to  segment  human-heart  cardiac  borders,  due  to  significant 
signal  loss,  noise,  and  non-uniformity  of  regional  intensities  in  ultra  sound  images. 
The  purpose  of  the  model  is  to  segment  two-chamber  end-systolic  endocardial 
human-heart  ultrasound  images  using  a prior  shape  and  a prior  intensity  profile. 

Feng,  Chen,  et  al.  developed  an  algorithm  to  generate  a prior  mean  shape 
by  clustering  extracted  from  eighty-five  two-chamber  end-systolic  endocardial 
human-heart  cardiac  borders.  These  eighty-five  contours  were  grouped  into  three 
clusters.  The  proposed  model  is  tested  against  the  mean  contour  from  the  cluster 
which  contained  thirteen  contours.  Prior  mean  shape  C * and  intensity  information 
I*  of  these  thirteen  contours  and  their  corresponding  images  which  have  0.62mm 
resolution  in  each  pixel  are  created  by  using  the  method  in  section  2.3.1  and  2.3.2. 
For  the  details,  refer  to  Chen,  Huang,  et  al.  [12].  The  prior  intensity  information 
and  a prior  shape  are  in  Figure  2-1.  The  left  image  in  Figure  2-1  represents  the 
average  image  as  prior  intensity  information.  The  solid  line  in  the  right  image  in 
Figure  2-1  is  the  average  curve  as  a prior  shape. 

To  show  the  effectiveness  of  the  local  non-rigid  deformation  term  in  the  model, 
a global  rigid  transformation  model  is  created  without  non-rigid  deformation  term 

in  Equation  (2-9).  Numerical  calculation  is  done  without  the  non-rigid  term. 

o validate  experimental  results  of  the  model,  distance  measurements  from  the 
results  of  the  model  with  non-rigid  term  from  expert’s  borders  from  Figure  2-4  to 
Figure  2-7  are  followed  in  the  Table  2-2.  Distance  measurements  to  the  results  of 
the  model  with  rigid  term  only  from  expert’s  borders  from  Figure  2-4  to  Figure  2-7 
are  presented  in  the  Table  2-3.  Tables  2-2  and  2-3  provide  the  effectiveness  of  the 
proposed  model. 
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In  Figures  (a)  2-2,  2-4,  2-5,  2-6,  and  2-7,  a solid  line  is  the  initial  contour 
in  a novel  image.  A prior  shape  is  used  as  an  initial  contour.  From  Figure  2-2 
to  Figure  2-3,  the  effectiveness  of  each  term  is  shown  in  the  model.  The  solid 
line  from  Figure  2-2  to  Figure  2-3  is  the  proposed  model  result  with  special 
condition  of  A*.  where  i = 1,2, 3, 4.  The  expert  endocardial  border  is  compared 
as  a dotted  line.  Figure  (b)  in  2-2  shows  the  model  performance  with  AL  = 0 and 
the  distance  measurement  from  expert’s  border  is  29. 1537(18. 7529mm).  Figure  (c) 
in  2-2  shows  the  model  performance  with  A2  = 0 and  the  distance  measurement 
from  expert’s  border  is  2.4027  (1.4897mm).  Figure  (a)  in  (2-3)  shows  the  model 
performance  with  A3  = 0 and  the  distance  measurement  from  expert  s border  is 
9.7669  (6.0555mm).  Figure  (b)  in  (2-3)  shows  the  model  performance  with  A4  = 0 
and  the  distance  measurement  from  expert’s  border  is  1.8930  (1.1737mm).  From 
Table  2-1,  it  is  easily  seen  that  the  major  contribution  terms  in  the  model  are 
the  first  and  third  terms.  From  Figure  2-4  to  Figure  2-7,  the  model  comparisons 
between  the  local  non-rigid  term  and  without  the  local  non-rigid  term.  For  Figure 
(b)  in  (2-4)  to  (2-7),  the  solid  line  is  the  model’s  result  with  the  local  non-rigid 
deformation  term.  The  expert  endocardial  border  is  represented  by  a dotted  line. 
For  Figure  (c)  in  2-4  to  2-7,  the  solid  line  is  the  result  of  the  model  without  a 
local  non-rigid  deformation  term  in  the  Equation  (2-9).  In  a similar  way,  the 
expert  endocardial  border  is  represented  as  a dotted  line.  From  the  Figures  and 
the  distance  measurement  results  from  the  expert’  border  in  Table  2-2  and  Table 
2-3.  it  is  easily  seen  that  the  performance  of  the  suggested  model  is  better  than 
the  model  without  a local  non-rigid  deformation  in  the  Equation  (2-9).  In  Figure 
2-4.  the  distance  result  from  the  model  with  non-rigid  term  is  1.8442(1. 1434mm) 
which  is  better  than  the  distance  result  2. 2160(1. 3739mm)  from  the  global  rigid 
transformation  model.  In  Figure  2-5,  the  distance  result  from  the  model  with 
non-rigid  term  is  1.8431(1. 1427mm)  which  is  better  than  the  distance  result 
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1.9670(1. 2195mm)  from  the  global  rigid  transformation  model.  In  Figure  2-6. 
the  distance  result  from  the  model  with  non-rigid  term  is  2. 2301(1. 3826mm.) 
which  is  better  than  the  distance  result  2. 8891(1. 7912mm)  from  the  global  rigid 
transformation  model.  In  Figure  2-7,  the  distance  result  from  the  model  with 
non-rigid  term  is  2. 0686(1. 2825mm.)  which  is  better  than  the  distance  result 
3. 0470(1. 8891mm)  from  the  global  rigid  transformation  model.  Since  the  local 
deformation  is  assumed  small  in  the  model,  the  distance  difference  between  the 
model  with  non-rigid  term  and  the  global  rigid  transformation  only  model  is  not 
big.  However,  Table  2-2  and  Table  2-3  shows  the  effectiveness  of  the  suggested 
model.  Finally,  Table  2-4  and  Table  2-5  presented  the  parameters  which  were 
chosen  as  examples  where  accurate  results  of  the  proposed  model  were  obtained. 

During  the  numerical  experiments,  the  parameters  Aj,  i = 1,2, 3, 4 are 
estimated  by  trial  and  error  basis.  There  are  several  different  ways  to  estimate 
unknown  parameters.  One  of  the  approaches  is  to  incorporate  unknown  parameters 
into  the  energy  function.  Then  find  steady  state  solutions  of  Euler-Lagrange 
equations  involved  with  these  parameters.  Statistical  approaches  have  been  also 
developed  to  get  best  estimation  of  unknown  parameters.  The  estimation  of 
parameters  can  be  presented  in  Gauss-Markoff  model  and  several  approaches 
including  least  squares  method  and  maximum-likelihood  method  have  been 
proposed  to  estimate  Gauss-Markoff  model.  For  the  details,  refer  to  [39]. 

A standard  technique  to  estimate  parameters  in  our  proposed  model  is 
introduced.  First,  fix  one  parameter,  say  A4  = 1.  If  we  have  a training  set  of  n 
contours  and  images  from  the  expert,  follow  the  procedure  from  section  2.3  to 
get  prior  shape  and  intensity  information.  Then  a,u,v,  /x,  R,T  can  be  measured 
using  the  discretized  solutions  of  Equations  (2-10)  through  (2-15).  Hence  n energy 
functions  Ej (a,,  Uj,  t»i, /ij,  Tj)  = 0,  i = 1,  2, ...,  n are  formed.  This  process 
results  in  the  following  matrix  equation  : 
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(a)  (b)  (c) 

Figure  2-2:  Subject  2 : Solid  Line  by  Our  Model  Result  and  Dotted  Line  by  Ex- 
pert (a)A  Initial  Contour,  (b)  Model  with  Ai  = 0,  (c)  Model  with  A2  = 0 


(a)  (b) 

Figure  2-3:  Subject  2 : Solid  Line  by  Our  Model  Result  and  Dotted  Line  by  Ex- 
pert (a)  Model  with  A3  = 0,  (b)  Model  with  A4  = 0 


An 

A21 

A„i 


A12 

A22 

A13 

A23 

Ai 

A2 

Bi 

b2 

An2 

An3 

1 

co 

<< 

1 

B„ 

(2-16) 


where  Aa  = \V{(G„)i*I*)(x,y)-a'V({Ga)i)*I){ui,vi)\2dx,  Ai2  = J)v(c.)(H2+ 

N 2)dx,  Ai 3 = /^.jdVuip+lVujpjcte,  and  B { = - fN(c.)  <5(dc*)9i(|V/|)(ui,  Vi)\Vdc> \dx, 


Denote  A := 

^11 

A21 

A12 

A22 

A13 

A 23 

, A := 

Ai 

A2 

, and  B := 

5i 

b2 

Ani 

An2 

An3 

A3 

Bn 

multiply  A 1 to  both  sides  of  Equation  (2-16).  Find  an  inverse  matrix  of  3 x 3 
matrix  A1  A and  multiply  to  both  sides  to  solve  the  equation  A* AX  = AlB.  Finally, 
desired  Ai,  A2,  and  A3  are  obtained. 
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Table  2-1:  Distance  Measurements  of  the  Result  of  the  Model  from  Expert’s  Bor- 
ders from  Figure  2-2  to  Figure  2-3 


Terms  in  Model 

(a)  Initial 

(b)  Result 

Subject  2 with  Ai  = 0 
Subject  2 with  A2  = 0 
Subject  2 with  A3  = 0 
Subject  2 with  A4  = 0 

29. 1273(18. 0589mm) 
29.1273(18.0589mm) 
29. 1273(18. 0589mm) 
29. 1273(18. 0589mm) 

29. 1537(18. 0753mm) 
2.4027  (1.4897mm) 
9.7669  (6.0555mm) 
1.8930  (1.1737mm) 

Figure  2-4:  Subject  1 : Solid  Line  by  Our  Model  Result  and  Dotted  Line  by  Ex- 
pert (a)An  Initial  Contour  (b)  Model  with  Non-Rigid  Term,  (c)  Model  with  Rigid 
Term  Only 


(a)  (b)  (c) 


Figure  2-5:  Subject  2 : Solid  Line  by  Our  Model  Result  and  Dotted  Line  by  Ex- 
pert (a)An  Initial  Contour  (b)  Model  with  Non-Rigid  Term,  (c)  Model  with  Rigid 
Term  Only 


Figure  2-6:  Subject  3 : Solid  Line  by  Our  Model  Result  and  Dotted  Line  by  Ex- 
pert (a)An  Initial  Contour  (b)  Model  with  Non-Rigid  Term,  (c)  Model  with  Rigid 
Term  Only 
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(a)  (b)  (c) 

Figure  2-7:  Subject  4 : Solid  Line  by  Our  Model  Result  and  Dotted  Line  by  Ex- 
pert (a)An  Initial  Contour  (b)  Model  with  Non-Rigid  Term,  (c)  Model  with  Rigid 
Term  Only 


Table  2-2:  Distance  Measurements  of  the  Result  of  the  Model  with  Non-Rigid 
Term  from  Expert’s  Borders  from  Figure  2-4  to  Figure  2-7 


Figures 

(a)Initial 

(b)Non-rigid  term  result 

Subject  1 
Subject  2 
Subject  3 
Subject  4 

25.1354  (15.5839mm) 
29. 1273(18. 0589mm) 
10.7297  (6.6524mm) 
4.4458  (2.7563mm) 

1.8442(1. 1434mm) 
1.8431  (1.1427mm) 
2.2301  (1.3826mm) 
2.0686  (1.2825mm) 

Table  2-3:  Distance  Measurements  of  the  Result  of  the  Model  with  Rigid  Term 
Only  from  Expert’s  Borders  from  Figure  2-4  to  Figure  2-7 


Figures 

(a)  Initial 

(b)Rigid  term  result 

Subject  1 
Subject  2 
Subject  3 
Subject  4 

25.1354  (15.5839mm) 
29. 1273(18. 0589mm) 
10.7297  (6.6524mm) 
4.4458  (2.7563mm) 

2.2160  (1.3739mm) 
1.9670  (1.2195mm) 
2.8891  (1.7912mm) 
3.0470  (1.8891mm) 

27 


Table  2-4:  The  Parameters  Used  for  Figure  2-4  through  Figure  2-7  with  Non-Rigid 
Term 


Figures 

(a)Al 

(b)A2 

(c)A3 

(d)A4 

Subject  1 

4.0 

0.0001 

500 

100 

Subject  2 

0.2 

0.0000001 

500 

100 

Subject  3 

0.2 

0.0000001 

500 

100 

Subject  4 

0.5 

0.0001 

500 

100 

Table  2-5:  The  Parameters  Used  for  Figure  2-4  through  Figure  2-7  with  Rigid 
Term  Only 


Figures 

(a)Al 

(b)A2 

(c)A3 

(d)A4 

Subject  1 

4.0 

0 

0 

100 

Subject  2 

0.2 

0 

0 

100 

Subject  3 

0.2 

0 

0 

100 

Subject  4 

0.5 

0 

0 

100 

2.8  Generating  Prior  Mean  Shape  and  Clustering 

Shape  analysis  has  become  one  of  the  most  important  and  interesting  areas 
and  has  a wide  range  of  applications  including  pattern  recognition,  object  recogni- 
tion, and  medical  image  analysis.  Due  to  more  advanced  algorithmic  models  that 
better  allow  researchers  to  analyze  images,  shape  analysis  is  now  a major  focus  to 
many  researchers  who  then  attempt  to  employ  mathematics,  statistics,  and  many 
other  computer  aided  techniques  to  further  theoretical  developments. 

Central  to  the  shape  analysis  problem  is  the  notion  of  a random  shape. 
Frechet  [23]  initiated  developments  of  this  area  for  analyzing  the  random  shapes, 
e.g.,  curves.  Matheron  [52]  and  David  Kendall  [34,  35,  36]  and  his  colleagues 
continued  after  him  for  mostly  further  theoretical  developments  of  the  analysis  of 
random  shapes. 

The  statistical  analysis  of  shapes  was  initiated  by  Sir  D'Arcy  Wentworth 
Thompson  [76].  A theory  and  practice  for  the  statistical  analysis  for  the  shapes 
have  been  developed  by  Bookstein  [4],  Dryden  and  Mardia  [20],  Carne  [6],  Kent 
and  Mardia  [37],  Cootes.  Taylor  and  colleagues  [18].  Their  ideas  are  tied  to  the 
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old  study  of  the  point-wise  representation  of  the  shapes:  objects  are  represented 
by  a finite  number  of  salient  points  or  landmarks  by  Frechet  [23],  Matheron  [52], 
and  David  Kendall  [34,  35,  36].  The  models  in  this  paper  follow  the  ideas  similarly 
from  Bookstein  [4],  Dryden  and  Mardia[20],  Carne  [6],  Kent  and  Mardia  [37],  and 
Cootes,  Taylor  and  colleagues  [18]  by  studying  the  point-wise  representation  of  the 
shapes  using  statistical  methods. 

One  of  the  objectives  of  the  shape  analysis  is  to  find  the  mean  shape  of 
the  given  contours  by  minimizing  a shape  related  energy  function.  Mean  shape 
plays  an  important  role  in  image  segmentation  and  registration  and  is  used 
as  prior  information  proposed  by  Chen,  leventon,  Paragios,  Yezzi,  Duncan,  et 
al.  [12,  13,  14,  15,  41,  42,  63,  64,  71,  85]  to  overcome  difficulties  which  include 
significant  signal  loss,  noise,  and  non-uniformity  of  regional  intensities.  To  find  the 
best  representation  of  the  given  contours,  several  statistical  methods  have  been 
developed  including  the  Procrustes  mean.  Principal  Component  Analysis  [73],  and 
generalized  orthogonal  Procrustes  analysis  [73].  The  suggested  model  is  not  only 
to  create  the  mean  shape  but  also  to  generate  alignments  between  the  mean  shape 
and  the  given  contours. 

The  other  goal  is  to  generate  the  mean  shape  and  clustering  of  the  given 
contours  using  the  Self-Organizing  map.  For  basic  pattern  recognition,  the  Self- 
Organizing  map,  which  was  developed  by  Professor  Teuvo  Kohonen  in  the  early 
1980s,  is  widely  used  and  can  be  visualized  as  a sheet-like  neural-network  array,  the 
cells  (or  nodes)  of  which  become  specifically  turned  to  various  input  signal  patterns 
or  classes  of  patterns  in  an  orderly  fashion.  The  learning  process  is  competitive  and 
unsupervised,  meaning  that  no  training  set  of  data  is  needed  to  define  the  correct 
output  (or  actually  the  cell  into  which  the  input  is  mapped)  for  an  input.  In  the 
basic  version,  only  one  map  node  (winner)  at  a time  is  activated  corresponding  to 
each  input.  The  details  of  the  self  organizing  map  can  be  found  in  [40]. 
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To  measure  the  similarities  between  the  cells(or  nodes),  the  Procrustes  method 
can  be  applied.  This  technique  was  given  a solid  mathematical  foundation  in 
the  seminal  work  of  Kendall  [35].  The  Procrustes  distance  is  a least-squares  type 
shape  metric  that  requires  one  to  one  correspondence  between  two  aligned  shapes. 
Several  Procrustes  methods  have  been  developed  including  an  area  difference 
distance  method  and  orthogonal  Procrustes  method.  The  orthogonal  Procrustes 
alignment  is  calculated  by  a point-wise  correspondence  and  there  is  no  need  for  the 
re-parametrization  process  to  obtain  the  correct  correspondence  between  contours. 
Meanwhile,  area  difference  measurement  method  needs  to  be  re-parametrized 
to  get  the  correct  correspondence  between  contours  in  each  iteration.  Hence 
the  orthogonal  Procrustes  method  has  an  advantage  of  faster  computation  time 
than  that  of  the  area  difference  distance.  In  this  dissertation,  these  two  methods 
combined  with  Self-Organizing  map  are  compared.  In  addition,  the  Self-  Organizing 
map  combined  with  Principal  Component  Analysis  is  presented  to  get  better  shape 
variations. 

2.8.1  Generating  Mean  Shape  with  Orthogonal  Procrustes  Measure- 
ment 

The  model  which  minimizes  a shape  related  energy  function  to  create  shape 
representation  with  the  orthogonal  Procrustes  method  is  presented.  This  model 
does  not  depend  upon  the  choice  of  the  initial  contour.  The  orthogonal  Procrustes 
analysis  finds  the  best  alignment  between  two  contours.  Since  the  method  is  fairly 
simple  and  calculation  time  is  fast,  it  is  widely  used.  The  followings  are  the  brief 
review  of  the  orthogonal  Procrustes  method  presented  by  Stegmann  and  Gomez 
[73]: 

2. 8. 1.1  Orthogonal  Procrustes  Method. 

The  Procrustes  distance  method  is  a least-square  type  shape  metric  that 
requires  two  aligned  contours  with  one-to-one  point  correspondence.  The  method 
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was  named  from  Procrustes(the  subduer),  son  of  Poseidon,  who  kept  an  inn 
benefiting  from  what  he  claimed  to  be  a wonderful  all-fitting  bed.  He  lopped  off 
excessive  limbage  from  tall  guests  and  either  flattened  short  guests  by  hammering 
or  stretched  them  by  racking.  The  victim  fitted  the  bed  perfectly  but,  regrettably, 
died.  To  exclude  the  embarrassment  of  an  initially  exact- fitting  guest,  variants  of 
the  legend  allow  Procrustes  two,  different-sized  beds.  Ultimately,  in  a crackdown 
on  robbers  and  monsters,  the  young  Theseus  fitted  Procrustes  to  his  own  bed. 

The  orthogonal  Procrustes  method  is  a least-square  shape  metric  that  requires  two 
aligned  shapes  with  one-to-one  correspondence.  The  alignment  part  involves  four 
steps  : 

STEP  1.  Compute  the  centroid  of  each  iteration. 

STEP  2.  Re-scale  each  shape  to  have  equal  size. 

STEP  3.  Align  with  respect  to  position  the  two  shapes  at  their  centroid. 

STEP  4.  Align  with  respect  to  orientation  by  rotation. 

2. 8. 1.2  Description  of  the  Suggested  Model 

A new  model  which  minimizes  a shape  related  energy  function  to  create 
mean  shape  with  the  orthogonal  Procrustes  method  is  presented  in  this  section. 
This  model  does  not  depend  upon  the  choice  of  the  initial  contour.  In  addition, 
the  model  generates  the  mean  shape  and  the  alignments  for  the  given  shapes  to 
this  mean  shape  simultaneously.  The  symbol  C*  is  the  mean  shape  of  the  given 
contours  C*i, . . . , Cn  and  Ri, . . . , Rn  are  the  rotation  matrices  with  respect  to 
the  angle  9 which  brings  the  best  fit  to  the  mean  shape  from  the  given  shapes, 
respectively.  Here  n is  the  total  number  of  given  contours.  The  model  finds 
C* . Ri, . . . , Rn  by  minimizing  the  energy  functional  : 

n 

E(C'.Rl,....Rn)  = J2\\c*  ~ R.QW2. 

i=l 


(2-17) 
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The  Equation  (2-17)  was  solved  by  finding  a steady  state  solution  of  the 
evolution  equations.  The  evolution  equations  are  associated  with  the  Euler- 
Lagrange  equations  of  the  Equation  (2-17).  Then  apply  the  gradient  descent 
method  to  the  evolution  equation  to  generate  the  mean  shape  C* . Then  use 
the  orthogonal  Procrustes  method  to  get  the  best  alignments  between  C*  and 
Ci, ... , Cn,  respectively.  This  gives  Ri,. . . ,Rn. 

The  iterative  approach  to  the  suggested  algorithm  is  the  following  six  steps  : 

STEP  1.  Compute  the  centroid  of  each  iteration. 

STEP  2.  Re-scale  each  shape  to  have  equal  size. 

STEP  3.  Choose  an  initial  estimate  of  the  mean  shape  (e.g.  the  arithmetic 
mean  of  the  given  n contours  ). 

STEP  4.  Align  all  the  remaining  shapes  to  the  mean  shape  to  get  R\, . . . , Rn.. 

STEP  5.  Recalculate  the  mean  shape  C*  by  applying  the  gradient  descent 
method  to  the  evolution  Equation  of  (2-17). 

STEP  6.  If  the  estimated  mean  has  changed,  return  to  STEP  4.  In  a similar 
way,  convergence  is  obtained  when  the  mean  shape  and  the  rotation  matrix  do  not 
change  within  an  iteration. 

2. 8. 1.3  Experimental  Results 

The  results  of  applications  to  the  eighty-five  human-heart  two-chamber  normal 
end-diastolic  epicardiac  borders  traced  by  experts  are  shown.  The  purpose  of 
the  model  is  to  generate  the  shape  representation  of  the  given  set  of  contours  by 
minimizing  a shape  related  energy  function.  The  model  does  not  depend  on  the 
choice  of  the  initial  contour  and  generates  the  mean  shape  and  alignments  to  the 
mean  shape  from  the  given  contours  simultaneously. 

Verification  of  the  results  are  done  by  three  different  measurements  : mean 
distance,  standard  deviation,  and  Tanimoto  measurement.  Let  = 1, 2,  • • • , n be 
given  contours  and  M be  an  average  shape  of  given  contours.  The  mean  distance 
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Figure  2-8:  Left  : A Initial  Contour,  Right  : The  Results  are  from  the  Shape  Re- 
lated Energy  Function  with  the  Mean  Shape  as  a Dotted  Line  with  85  human-heart 
end-diastolic  epicardial  borders 


is  L)i_i  d(AltM)  wkere  (i  js  a ^stance  between  A and  M . The  standard  deviation 

n ’ 1 

is  LAi  std(At,M) , wjiere  std(A,  M)  is  a standard  deviation  between  At  and  M. 
Finally,  the  Tanimoto  distance  is  ^=iTi^Ai'M^ , where  Td  — 1 


A,-M 


IIAIP+IIMP-Ai-M' 

This  measurement  is  used  to  compare  the  dissimilarity  of  the  given  contours  to  an 
average  shape.  The  details  of  the  Tanimoto  similarity  measurement  are  presented 
by  Kohonen  [40]. 

Figure  2-8  represents  the  shape  related  energy  function  model  result  applied 
eighty-five  human-heart  end-diastolic  epicardial  borders.  For  Figure  2-8,  a dotted 
line  in  the  left  one  is  the  initial  contour  of  the  given  contours.  The  dotted  line  in 
the  right  one  is  mean  shape  of  the  given  contours.  Three  distance  measurements, 
namely  mean  distance,  standard  deviation,  and  Tanimoto  measurement  of  the 
result  of  the  models  from  Figure  2-8  is  in  Table  2-6. 

Table  2-6:  Distance  Measurements  of  the  Result  of  the  Models  from  Figure  2-8, 
n=85 


Figure 

(a)Mean  Distance 

(b)  Standard  Deviation 

(b)Tanimoto  Distance 

Figure  2-20 

0.0635 

0.0535 

0.0088 
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2.8.2  Generating  Mean  Shape  and  Clustering  using  Self-Organizing 
Map 

2.8.2. 1 Description  of  the  Self-Organizing  Map(SOM) 

In  this  section,  the  Self-Organizing  map  is  described.  The  following  is  the  brief 
description  of  the  Self-Organizing  map  from  Kohonen  and  Honkela  [30,  40]: 

The  set  of  input  samples  is  described  by  a real  vector  x(f)  € Rn,  where  t 
is  the  index  of  the  sample,  or  the  discrete-time  coordinate.  Decide  number  of 
clusters  in  advance,  say  k.  Each  clustering  i — 1, 2, . . . , k in  the  map  contains  a 
model  vector  m,-(t)  6 Rn,  which  has  the  same  number  of  elements  as  the  input 
vector  x(t).  The  stochastic  Self-Organizing  Map  algorithm  performs  a regression 
process.  Thereby,  the  initial  values  of  the  components  of  the  model  vector,  mj(t), 
may  even  be  selected  at  random.  Any  input  item  is  thought  to  be  mapped  into 
the  location,  the  mi;(t)  of  which  matches  best  with  x(f)  in  some  metric.  The  self- 
organizing algorithm  creates  the  ordered  mapping  as  a repetition  of  the  following 
basic  tasks:  First,  compare  an  input  vector  x(t)  with  all  the  model  vectors  m i(t). 
Best-matching  unit  (node)  on  the  map,  i.e.,  the  node  where  the  model  vector  is 
most  similar  to  the  input  vector  in  some  metric  (e.g.  Euclidean)  is  identified.  This 
best  matching  unit  is  often  called  the  winner.  Then,  the  model  vectors  of  the 
winner  and  a number  of  its  neighboring  nodes  in  the  array  are  changed  towards 
the  input  vector  according  to  the  learning  principle  specified  below.  Adaptation  of 
the  model  vectors  in  the  learning  process  may  take  place  according  to  the  following 
equations: 

uii(t  + 1)  = m,(f)  + a(f)[x(t)  - m i(t)],if  i € Nc(t ) 

(2-18) 

m,(f+  1)  = nii(i),  otherwise, 

where  t is  the  discrete-time  index  of  the  variables,  the  factor  a(t)  € [0, 1]  is  a 
scalar  that  defines  the  relative  size  of  the  learning  step,  and  Nc(t)  specifies  the 
neighborhood  around  the  winner  in  the  map  array.  At  the  beginning  of  the  learning 
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process  the  radius  of  the  neighborhood  is  fairly  large,  but  it  is  made  to  shrink 
during  learning.  This  ensures  that  the  global  minimum  is  obtained  already  at 
the  beginning,  whereas  towards  the  end,  as  the  radius  gets  smaller,  the  local 
corrections  of  the  model  vectors  in  the  map  will  be  more  specific.  The  factor  a(i) 
also  decreases  during  learning. 

2. 8. 2. 2 Self-Organizing  Map  with  Procrustes  Methods 

The  purpose  of  the  model  is  to  generate  the  shape  representation  and  clas- 
sification of  the  given  contours  using  the  Self-Organizing  map  combined  with  the 
Procrustes  methods.  This  technique  has  the  advantage  of  generating  shape  rep- 
resentation of  each  cluster  and  classifying  given  contours  simultaneously.  In  this 
section,  the  suggested  model  is  explained.  The  following  is  the  description  of  the 
proposed  model: 

STEP  1.  First,  centralize  and  normalize  all  the  n training  contours  C*.  to 
get  the  equal  length.  If  the  n training  contours  Cj,  i = 1, ...  ,n  to  be  grouped 
into  k clusters,  take  k arbitrary  contours  as  the  initial  contours,  denoted  by  rrij{ 0) 

{j  = 

STEP  2.  At  ( t + 1)  iteration,  randomly  select  a contour  denoted  by  X(t  -t-  1) 
from  the  training  set,  and  compare  the  disparity  in  shape  between  X(t  - (-  1)  and 
each  of  rrij(t)  (j  = 1, . . . , k). 

STEP  3.  To  do  this  comparison,  align  X(t  + 1)  to  each  rrij(t),  and  denote  the 
aligned  X(t  + 1)  by  Xj(t  + 1)  = RjX(t  + 1)  + Tj. 

STEP  4.  Then  we  compute  A3  AD(Xj(t  + 1 ):rrij)  by  using  the  orthogonal 

Procrustes  method  in  (2. 8. 1.1). 

STEP  5.  Suppose  A\  is  the  smallest  number  in  A3 , j = 1,  ...,k.  Keep  m2(t)  , 

. . . , mfc(f)  unchanged,  and  update  mi(t)  by 

m\{t  + 1)  = mx(t)  + a(t)[Xi(f  + 1)  - mi(f)], 
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where  a(t)  is  a smooth  function  of  t , and  decreases  to  zero  as  t — » oo.  Then 
normalize  mi(t). 

STEP  6.  After  a large  number,  say  t = N,  k average  shapes  nij(N)  ( j = 

1, . . . , k)  are  generated. 

STEP  7.  Then  k clusters  are  formed  by  the  contours  that  are  closest  to  the 
average  shapes.  The  closeness  is  again  measured  by  an  orthogonal  Procrustes 
measurement  or  area  difference  distance  measurement. 

2. 8. 2. 3 Self-Organizing  Map  Combined  with  Principal  Component 
Analysis 

The  purpose  of  this  model  is  to  generate  the  mean  shape  and  clustering  of 
the  given  contours  using  the  Self-Organizing  map  combined  with  the  Principal 
Component  Analysis  method.  The  method  is  very  similar  to  Section  (2. 8. 2. 2).  To 
get  better  results,  the  given  contours  are  pre-aligned  using  the  method  in  Section 
(2. 8. 1.1).  Then  the  algorithm  from  Section  (2. 8. 2. 2)  is  applied.  After  that  Principal 
Component  Analysis  is  applied  to  each  group  of  clusters.  The  following  is  the  basic 
summary  of  the  proposed  model: 

STEP  1.  Pre-align  the  given  contours  by  applying  the  orthogonal  Procrustes 
algorithm  from  Section  (2. 8. 1.1). 

STEP  2.  Apply  the  STEP  1.  through  STEP  7.  from  Section  (2. 8. 2. 2)  to  get 
the  mean  shape  and  clustering  to  the  pre-aligned  contours. 

STEP  3.  Update  the  mean  shape  and  clustering  of  each  group  by  applying  the 
Principal  Component  Analysis. 

2. 8. 2. 4 Experimental  Results 

The  results  of  applications  to  the  eighty-five  human-heart  two-chamber  normal 
end-diastolic  epicardiac  borders  are  shown.  The  model  does  not  depend  upon  the 
choice  of  the  initial  contour  and  generates  the  mean  shape  and  clustering  from  the 
given  contours  simultaneously. 
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Verification  of  the  results  are  done  by  three  different  measurements;  mean 
distance,  standard  deviation,  and  Tanimoto  measurement.  Let  Ai,  i = 1,  2,  • • ■ , n be 
given  contours  and  M be  an  average  shape  of  given  contours.  The  mean  distance 
js  ^ where  d is  a distance  between  and  M . The  standard  deviation 
js  E,=i  std(A,,M) , w}iere  std(At.M)  is  a standard  deviation  between  A,  and  M. 
Finally,  the  Tanimoto  distance  is  E.‘=i  TfPVA/)  where  Td  = 1 - V ¥ , ,, 

This  measurement  is  used  to  compare  the  dissimilarity  of  the  given  contours  to  an 
average  shape.  The  detail  of  the  Tanimoto  similarity  measurement  is  in  [40]. 

From  Figure  2-9  to  Figure  2-13  are  the  results  which  were  obtained  when 
this  method  was  applied  to  the  eightv-five  human-heart  two-chamber  normal 
end-diastolic  epicardiac  borders  to  obtain  clustering  and  the  mean  shape  of  each 
cluster.  The  self-organizing  map  is  used.  The  closeness  was  measured  by  the  area 
difference  method  or  the  orthogonal  procrustes  method.  In  the  experiment,  we  use 
o(f)  = 0.9(1  — t/N),  where  t is  the  time  step  and  N is  total  number  of  iterations. 
The  function  a is  smooth  and  decreases  to  zero  as  t approaches  N.  Figure  2-9  is 
the  numerical  results  using  Self-Organizing  map  with  the  area  difference  distance 
method.  In  a similar  way,  Figure  2-10  is  the  numerical  results  using  Self- Organizing 
map  with  orthogonal  Procrustes  method.  Since  orthogonal  Procrustes  method  is 
using  one-to-one  correspondence,  we  do  not  need  the  re-parameterizing  process 
to  find  the  correct  correspondence  which  make  the  calculation  time  faster  in  each 
step  than  that  of  using  area  difference  distance  measurement.  It  is  easily  seen  that 
the  orthogonal  method  is  more  efficient  than  area  difference  method  from  Table 
2-7.  To  get  better  mean  shape,  Figure  2-11  shows  the  numerical  results  using  the 
method  fiom  section  2. 8. 2. 3.  Pre-alignments  and  applying  principal  component 
analysis  give  better  results  which  can  be  easily  seen  in  Table  2-7.  From  Figure  2-12 
and  Figure  2-13,  numerical  results  with  an  application  to  the  eighty-five  human- 
heart  two-chamber  normal  end-diastolic  epicardiac  borders  to  cluster  into  three 
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Figure  2-9:  Two  Chamber  ED  EPI  Normal  85  Patients  Heart  Contours  with  5000 
Iterations  using  SOM  + Area  Difference  Method.  The  Dotted  Line  Represents  the 
Average  of  the  Cluster  with  n=85 


Figure  2-10:  Two  Chamber  ED  EPI  Normal  85  Patients  Heart  Contours  with  5000 
Iterations  using  SOM  + orthogonal  Procrustes  method.  The  Dotted  Line  Repre- 
sents the  Average  of  the  Cluster  with  n=85 


groups.  Figure  2-12  shows  the  results  using  the  method  from  section  2. 8. 2. 2 with 
orthogonal  Procrustes  method.  First  group,  there  are  30  contours,  28  contours  are 
in  the  second  group.  Finally  27  contours  are  in  the  third  group.  Numerical  results 
using  the  method  from  section  2. 8. 2. 3 is  in  Figure  2-13.  Distance  measurements 
comparisons  from  Table  2-8  and  Table  2-9  show  the  effectiveness  of  the  suggested 
models. 
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Figure  2-11:  Two  Chamber  ED  EPI  normal  85  Patients  Heart  Contours  with  5000 
Iterations  using  SOM  + PCA.  The  Dotted  Line  Represents  the  Average  of  the 
Cluster  with  n=85 

Table  2-7:  Distance  Measurements  of  the  Result  of  the  Models  from  Figure  2-9  to 
Figure  2-11.  Md  Represents  Mean  Distance.  STDd  Represents  Standard  Deviation. 
Td  Represents  Tanimoto  Distance  with  n=85 


Figures 

(a  )Md 

(b  )STDd 

(c  )Td 

Figure  2-21 

0.4915 

0.0504 

0.7036 

Figure  2-22 

0.0063 

0.0054 

0.0082 

Figure  2-23 

0.0057 

0.0051 

0.0044 

Figure  2-12:  Two  Chamber  ED  EPI  Normal  85  Patients  Heart  Contours  with  5000 
Iterations  into  3 Groups  using  SOM+Orthogonal  Procrustes  Measurement.  The 
Dotted  Line  Represents  the  Average  of  the  Cluster  with  n=85.  Left  : 30  contours, 
Middle  : 28  contours,  Right  : 27  contours 

Table  2-8:  Distance  Measurements  of  the  Result  of  the  Models  from  Figure  2-12. 
Md  Represents  Mean  Distance.  STDd  Represents  Standard  Deviation.  Td  Repre- 
sents Tanimoto  Distance 


Groups 

(a  )Md 

(b  )STDd 

(c  )Td 

Group  l(ni  = 30) 

0.0409 

0.0128 

0.0605 

Group  2 (n2  = 28) 

0.0426 

0.0128 

0.0529 

Group  3(n3  = 27) 

0.0286 

0.0103 

0.0350 
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Figure  2-13:  Two  Chamber  ED  EPI  Normal  85  Patients  Heart  Contours  with 
5000  Iterations  using  SOM  + PCA.  The  Dotted  Line  represents  the  Average  of  the 
Cluster  with  n=85.  Left  : 30  contours,  Middle  : 28  contours,  Right  : 27  contours 


Table  2-9:  Distance  Measurements  of  the  Result  of  the  Models  from  Figure  2-13. 
Md  Represents  Mean  Distance.  ST Dd  Represents  Standard  Deviation.  Td  Repre- 
sents Tanimoto  Distance 


Groups 

(a  )Md 

(b  )STDd 

(cm 

Group  l(ni  = 30) 

0.0345 

0.0098 

0.0520 

Group  2(712  = 28) 

0.0415 

0.0101 

0.0505 

Group  3(n3  = 27) 

0.0265 

0.0097 

0.0306 

CHAPTER  3 

A REGION-BASED  MODEL  FOR  IMAGE  SEGMENTATION  AND 

REGISTRATION 

3.1  Image  Segmentation  Using  a Modified  Mumford-Shah  Technique 
3.1.1  Introduction 

Capturing  the  boundary  of  the  given  object  is  of  primary  importance  in  image 
processing,  with  applications  to  computer  vision  and  medical  imaging.  Region 
based  segmentation  techniques  have  been  developed  by  using  the  region  statistics 
or  homogeneity  proposed  by  Chan,  Vese,  Cremers,  Tsai,  Mumford,  Shah,  Yezzi, 

Zhu,  et  al.  [9,  10,  19,  22,  57,  75,  87,  93], 

The  most  celebrated  region  based  image  segmentation  model  is  introduced  by 
Mumford  and  Shah  in  1989  [57].  Using  this  model  for  a segmentation,  an  image  is 
decomposed  into  a set  of  regions  within  the  bounded  open  set  Q and  these  regions 
are  separated  by  smooth  edges  T.  Let  Q C Rn  be  a bounded  domain  with  Lipschitz 
boundary  and  u : Cl  — > R.  The  model  is  formulated  as  a variational  problem  to 
minimize  the  following  functional  : 

min  MS(u,  T)  = A f (u  - ifdx  + fiH(T)  + [ \Vu\2dx,  (3-1) 

“T  Jn  Jn\ r 

where  //stands  for  (n  — 1)  dimensional  Hausdorff  measure  and  A and  fi  are 

parameters  which  balances  three  terms  in  the  model.  Here  T is  to  be  a closed 

subset  of  Q given  by  the  union  of  a finite  number  of  curves.  It  represents  the  set 

of  ’’ edges”  (i.e.  boundaries  of  homogeneous  regions)  in  the  given  image  I.  The 

function  u is  the  piecewise  smooth  approximation  to  I.  Because  of  the  last  term  in 

Equation  (3-1)  taken  over  the  set  fi/T,  u is  smooth  in  each  connected  component  of 

n/r. 
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A simplification  of  Equation  (3-1)  is  to  find  contour  which  separate  two 
homogeneous  regions  with  average  intensity  c\  and  c2.  This  model  is  referred  to 
as  the  piecewise  constant  Mumford-Shah  model.  Hence  u = Clip  + c2ln/r;  where 
T C Q.  In  this  case,  Equation  (3-1)  reduces  to  the  following  minimization  problem  : 

min  £'(c1,c2,r)  = A / (u  - c{)2dx  + A / (u  - c2)2 dx  + fj,Length(T) , (3-2) 

Cl’C2’r  Jr  Jn/r 

where  A and  p are  parameters  which  balances  three  terms  in  the  model.  Model 

(3-2)  is  known  as  the  two-phase,  piecewise  constant  Mumford-Shah  model.  Chan 

and  Vese  [9]  reformulate  this  model  in  the  context  of  region  based  active  contour  : 

min  E(ci,c2,r)  = A / (it  - Cifdx  + A / (it  - c2)2dx  + uLengthiT),  (3-3) 
Cl-C2-r  J n+  Jn_ 

where  T partitions  Q into  interior  fl+  and  and  A and  p are  parameters  which 
balances  three  terms  in  the  model.  Chan  and  Vese  proposed  a piecewise  constant 
Mumford-Shah  model  in  [9,  10]  by  using  a level  set  formulation  presented  by  Osher 
and  Fedkiw  [58].  In  this  approach,  the  boundary  T is  represented  as  the  zero  level 
set  of  a Lipschitz  continuous  function  0 : f 1 -*  i?.  The  purpose  of  this  model  is  to 
express  Equation  (3-3)  in  terms  of  the  level  set  function  0 : 


min  E(a,c2,r)  = A f H(<p)(u-cx)2dx  + A f (1  - H((f)))(u-c2)2dx+  f \VH(<p)\dx, 
C1,C2’  Jn  Jn  J n 

(3-4) 

where  A is  a parameter  which  balances  three  terms  in  the  model  and  H(<f>)  : R — > R 
is  the  Heaviside  function  such  that 


H(ct>) 


0,  if  0 < 0, 

1,  0. 


Multiphase  frameworks  have  also  been  developed  by  Chan  and  Vese  [79].  Devel- 
opments of  variational  level  set  implementation  techniques  based  Mumford-Shah 
model  are  followed  by  Esedoglu,  Gibou,  Likar,  et  al.  [22,  26,  43]. 
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In  the  theory  of  T-convergence  [1,  21,  49],  the  measuring  an  edge  T length 
term  in  the  Mumford-Shah  model  can  be  approximated  by  a quadratic  integral  of 
an  edge  signature  function  p(x)  such  that 

/ (e|  Vp|2  + —)dx,  e <C  1 

J n 4e 

by  Ambrosio  and  Tortorelli  in  1990  [1].  The  segmentation  is  represented  by 
characteristic  function  using  phase  fields  in  Esedoglu,  Likar,  et  al.  [22,  43].  A 
diffuse  interface  approximation  for  Equation  (3-2)  proposed  by  Esedoglu,  et  al.  [22] 
is  given  by  the  following  sequence  of  energies  : 

MSe(u,  Ci,  c2)  = [ e|Vu|2  + -W(u)  + \{u2(I  - cq)2  + (1  - u)2(I  - c2)2}dx,  (3-5) 
Jn  £ 

where  e > 0,  A is  a parameter,  u is  the  characteristic  function,  and  the  potential 
W(u)  — u2(l  — u)2.  This  potential  function  improves  sharpness  of  edges.  The 
details  of  phase  field  theory  can  be  found  in  Baldo,  Modica,  Shen,  Zhou,  et  al. 

[3,  55,  70,  82], 

Recently,  fuzzy  segmentation  algorithm  using  a modified  Mumford-Shah  model 
is  proposed  by  Shen  [70]  and  Thiruvenkadam  [77],  Shen  [70]  and  Thiruvenkadam 
[77]  use  0 < p < 1 function  instead  of  a characteristic  function  u in  Equation  (3-5). 
Shen  [70]  proposed  the  following  minimization  problem  : 


min  E(p,ci,c2)  = X / p{I  - cx)2  + (1  - p)(I  - c2fdx+ 

P,C1,C2  JQ 

! (9e(|Vp|2  + | V(1  - p)|2)  + 2(p(1~P))2),tt.  (3-6) 

Jn  e 

where  e > 0 and  p is  a probability  measure  of  image  pattern  with  0 < p < 1.  In 
[77],  p = cos2 {9)  is  presented,  where  0 < 9 < | : 


min  E(6,  c\,c2)  = [ cos2(d)(I 

0,C1,C2  Jn 


- Cl)2  + sin2{6){I  - c2)2  + A|V6>|2  + ^(|V//(<l>)|)di, 


(3-7) 
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where  0 = cos2(6)  — sin2(6)  and  A and  )i  are  parameters. 


In  this  section,  a region  based  variational  partial  differential  equation  model 
for  an  image  segmentation  is  presented  with  an  application  to  synthetic  and 
simulated  human  brain  MRI  images.  The  purpose  of  the  model  is  to  segment 
borders  of  the  given  images.  The  segmentation  is  obtained  by  using  a modified 
Mumford-Shah  segmentation  technique.  The  proposed  algorithm  is  motivated 
by  Esedoglu  [22],  Shen  [70],  Thiruvenkadam  [77],  et  al..  A fuzzy  segmentation 
in  our  model  is  similar  to  Shen  [70]  and  Thiruvenkadam  [77],  but  the  length  of 
boundary  term  of  proposed  model  is  different  from  Thiruvenkadam  [77].  The 
length  of  boundary  term  is  motivated  by  Esedoglu  [22],  but  our  model  is  variational 
which  differs  from  Esedoglu  [22],  Our  model  is  applied  to  the  synthetic  image  and 
simulated  noisy  human  brain  images.  Numerical  experiments  show  the  effectiveness 
of  the  proposed  model  to  segment  given  images  and  the  robustness  to  a noise.  This 
section  is  organized  as  follows:  In  section  3.1.2,  a modified  Mumford-Shah  model  is 
described.  Euler-Lagrange  equations  of  the  suggested  model  are  presented  in  this 
section  as  well.  In  section  3.1.3,  numerical  methods  pertaining  the  proposed  model 
are  explained.  Experimental  results  of  the  model  which  were  applied  to  synthetic 
data  and  simulated  human  brain  MRI  images  are  shown  in  section  3.1.4. 

3.1.2  Description  of  a Modified  Mumford-Shah  Model 

In  this  section,  a region  based  segmentation  is  introduced.  The  segmentation 
is  attained  by  using  a modified  Mumford-Shah  segmentation  technique.  The  two 
phase  case  is  assumed  for  the  simplicity  of  the  problem  in  our  model.  The  model  is 
aimed  to  find  9,  C\ , and  C2  by  minimizing  the  energy  functional: 


(3-8) 
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where  I is  a novel  image,  fi  is  domain,  and  A and  p are  parameters  balancing 
the  influences  from  the  two  terms  in  the  model.  Here,  0 < 9 < n/2. 

The  first  term  forces  sin2(9(x)),  towards  1 if  I(x ) is  close  to  C\  and  towards  0 
if  I(x)  is  different  from  ci,  for  every  x E fb  In  a similar  way,  cos2  (9(x)) , towards 
1 if  I(x)  is  close  to  c2  and  towards  0 if  I(x)  is  different  from  C2,  for  every  x € Cl. 
The  second  term  is  the  length  of  boundary  term.  For  every  x € Cl,  denote 
<f)(x)  = cos2(9(x))  — sin2(9(x)).  If  (j>  > 0,  then  assign  1 to  inside  of  the  given  shape 
c.  If  < 0,  then  assign  0 to  inside  of  the  given  shape  c.  Finally,  the  segmented 
binary  image  is  formed. 

3. 1.2.1  Euler-Lagrange  Equations  of  the  Model 

The  evolution  equations  associated  with  the  Euler-Lagrange  equations  in 
Equation  (3-8)  are 


— = — \{sin(29(x))(I(x)  — c{)2  — sin(29(x))(I(x)  — c2)2}— 

{2esm(4#(ir))|V0|2  — 2esin2(29(x))A9+ 
fj,isin39(x)cos39cos(29(x)) . 

ZTl 

d6 

— = 0,  on  dCl 
on 


(3-9) 


= —A  f sin29(x)(I(x)  — c{)dx , (3-10) 

ot  Jn 

= — A f cos29(x)(I(x)  — c2)dx.  (3-11) 

ut  Jn 

3.1.3  Numerical  Methods 

In  this  section,  numerical  methods  to  solve  the  Equations  (3-9)-(3-ll)  are 
explained.  (3-8)  was  solved  by  finding  a steady  state  solution  of  (3-9)-(3-ll).  To 
discretize  the  equations  (3-9)-(3-ll),  we  use  a finite  difference  scheme  and  we 
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applied  the  gradient  descent  method.  Let  denote  the  value  of  9 at  the  pixel  X;  = 
ih  and  yj  = jk,  where  h and  k is  the  pixel  size.  Denote  9(xi,yj,tn)  by  9”..  The  time 
derivative  9t  at  (i,  j,tn)  is  approximated  by  the  forward  difference  scheme: 

/an+l  an 

et(ij,tn)=  ij  a i 


where  At  is  the  time  step.  Moreover,  (3-2)-(3-4)  were  implemented  by  the  following 
iteration  scheme  : 


(*.)?,  = 


fin  an 


an  fin 


(*»)«  = 
(«  = 


/in  _i_  fin 

h2 

k2 


3.1.4  Experimental  Results 

The  goal  of  this  experiment  is  to  segment  the  boundary  of  the  object  in  image 
using  a modified  Mumford-Shah  technique.  Throughout  the  experiments,  A = 1 
and  /i  = ~ are  used.  The  numerical  results  with  an  application  to  the  synthetic 
image  are  from  Figure  3-1  to  Figure  3-3.  Left  image  of  Figure  3-1  represents  a 
novel  image,  middle  image  is  9(x),  and  the  right  image  is  H(4>(x))  = 1,  if  4>(x)  > 0 
and  H(<f>(x))  = 0,  if  </>(x)  < 0,  where  <t>(x)  — cos29(x)  — sin29(x)  and  for  every 
I 6 fi.  In  Figure  3-2,  left  one  is  the  given  image  and  right  one  is  the  boundary  of 
the  segmented  image.  Figure  3-3  represents  values  of  c\  and  C2  in  every  iteration. 
Initial  value  of  C\  is  0.9  and  converges  to  0.9926  after  80  iterations.  In  a similar 
way,  C2  converges  to  0.1925  from  an  initial  value  0.5  after  80  iterations. 

The  results  of  an  application  to  the  simulated  human  brain  MRI  images  are 
from  Figure  3-4  to  Figure  3-9.  The  data  is  obtained  from  http://www.bic.mni. 
mcgill . ca/brainweb.  Simulated  MRI  volumes  for  normal  brain  with  T1  modality 


Figure  3-1:  Left  : A Novel  Image,  Middle  : 0(x),x  € Q,  Right  : Segmented  Image 
by  H(c/)(x)),x  € 0 


Figure  3-2:  Left  : A Novel  Image,  Right  : The  Boundary  of  the  Segmented  Image 


10  30  X 40  SO 


Figure  3-3:  Left  : C\  values,  Right  : c2  values  in  each  iteration  from  Figure  3-1 
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Figure  3-4:  Left  : A Novel  Image,  Middle  : 8(x),x  € O,  Right  : Segmented  Image 
by  H((f)(x)),x  € f i 


Figure  3-5:  Left  : A Novel  Image,  Right  : The  Boundary  of  the  Segmented  Image 

are  obtained.  Slice  thickness  is  1 mm  with  9%  noise  and  20%  non-uniformity 
intensity  level.  Numerical  results  from  our  model  show  From  Figure  3-4  to  Figure 
3-6.  The  goal  is  to  segment  light  blue  color  of  the  given  simulated  brain  image. 

Left  image  of  Figure  3-4  represents  9{x)  values  and  the  right  image  is  H{<fi(x))  = 1, 
if  (f>(x)  > 0 and  H((f>(x ))  = 0,  if  (j)(x)  < 0,  where  (f>{x ) = cos29{x)  — sin29(x)  and  for 
every  x G Cl.  The  left  one  is  given  image  and  the  right  one  is  the  boundary  of  the 
segmented  image  in  Figure  3-5.  Figure  3-6  represents  values  of  Ci  and  c-i  in  every 
iteration.  Initial  value  of  C\  is  0.01  which  converges  to  0.2126  after  80  iterations  and 
c2  converges  to  0.6467  from  an  initial  value  0.4  after  80  iterations. 

Another  numerical  result  with  an  applications  to  simulated  human  brain  image 
are  shown  from  Figure  3-7  to  Figure  3-9.  The  goal  is  to  segment  light  yellow  color 
of  the  simulated  brain  image.  Left  image  of  Figure  3-7  represents  6{x)  values  and 
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Figure  3-6:  Left  : C\  values,  Right  : c2  values  in  each  iteration  from  Figure  3-4 


Figure  3-7:  Left  : A Novel  Image,  Middle  : 6(x),x  6 Q,  Right  : Segmented  Image 
by  H(<j)(x)),x  € fl 

the  right  image  is  H(cj)(x))  = 1,  if  <f>(x)  > 0 and  H(<j)(x ))  = 0,  if  cj)(x)  < 0,  where 
(f)[x)  — cos26(x)  — sin29(x)  and  for  every  x G Q.  The  left  one  is  given  image  and 
the  right  one  is  the  boundary  of  the  segmented  image  in  Figure  3-8.  Figure  3-9 
represents  values  of  C\  and  c2  in  every  iteration.  Initial  value  of  Ci  is  0.01  which 
converges  to  0.1110  after  225  iterations  and  c2  converges  to  0.3474  from  an  initial 
value  0.4  after  225  iterations. 

3.2  Region  Based  Segmentation  with  Shape  Prior 
3.2.1  Introduction 

There  have  been  several  algorithms  to  develop  region  based  segmentation 
technique  using  the  region  statistics  or  homogeneity  presented  by  Chan,  Vese, 
Mumford,  Shah,  Zhu,  et  al.  [9,  10,  57,  93] 
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Figure  3-8:  Left  : A Novel  Image,  Right  : The  Boundary  of  the  Segmented  Image 


Figure  3-9:  Left  : c\  values,  Right  : ci  values  in  each  iteration  from  Figure  3-7 

Due  to  significant  signal  loss,  noise,  and  non-uniformity  of  regional  intensities, 
prior  information  is  necessary  to  make  the  resulting  segmentation  more  accurate. 
Cremers,  et.  al.  [19]  incorporated  statistical  shape  knowledge  into  Mumford- 
Shah  segmentation  scheme  [57]  by  minimizing  a functional  that  includes  both 
shape  energy  and  the  Mumford-Shah  energy.  A parametric  shape  model  was 
built  proposed  by  Tsai,  et  al.  [75]  and  the  parameters  adjusted  to  minimize 
a region  based  objective  function  which  provides  the  segmentation.  In  Chen, 
Tagare,  et  al.  [15],  the  shape  model  was  created  by  averaging  the  aligned  training 
samples.  Better  maps  were  gathered  by  matching  the  evolving  interface  to  the 
shape  prior.  The  geometric  active  contour  was  minimized  and  the  shape  disparity 
energy  function  between  the  active  contour  and  prior  shape  was  measured  by  a 
distance  function.  Recently,  deformable  models  incorporated  by  the  prior  shape 
and  intensity  information  have  been  proposed  by  Metaxas,  Huang,  et  al.  [31,  32], 
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A class  of  deformable  models  called  MetaMorphs  which  possess  both  shape  and 
interior  texture  statistics  integrate  boundary  and  region  information  coherently  in 
a common  variational  framework  [32],  The  MetaMorphs  models  [32]  are  extended 
into  a unified  shape  and  intensity  space  by  implicitly  representing  shapes  as  images 
in  the  space  of  distance  transforms.  A novel  stochastic  chord  based  matching 
algorithms  align  global  transformation  considering  both  shape  and  gray-level 
intensity  information  [31].  In  the  proposed  model,  it  is  assumed  that  the  prior 
shape  is  close  to  a novel  image  which  is  a similar  approach  from  An,  Chen,  et  al. 

[2].  But  the  presented  model  is  designed  for  the  images  which  have  homogenous 
regions  and  it  is  simple  to  implement. 

In  this  section,  a new  variational  image  segmentation  model  based  on  the 
region  intensity  information  using  shape  knowledge  is  suggested.  The  goal  of  the 
model  is  to  develop  the  shape  knowledge  based  segmentation  technique.  The  prior 
shape  information  is  obtained  by  using  a modified  Mumford-Shah  segmentation. 
The  prior  shape  information  is  supporting  the  image  segmentation  process  to  a 
given  novel  image.  But  prior  shape  information  is  formed  a binary  image  in  our 
model.  This  section  is  organized  as  follows:  In  section  3.2.2,  the  proposed  model 
is  introduced.  Euler-Lagrange  equations  of  the  suggested  model  are  presented  in 
this  section  as  well.  In  section  3.2.3,  numerical  methods  pertaining  the  proposed 
model  are  explained.  Experimental  results  of  the  model  which  were  applied  to  2D 
synthetic  data  and  human  brain  MRI  images  are  shown  in  section  3.2.4. 

3.2.2  Description  of  the  Proposed  Model 

A new  variational  image  registration  model  based  on  the  region  intensity 
information  using  shape  knowledge  is  suggested.  The  goal  of  the  model  is  to 
develop  the  prior  shape  knowledge  based  image  segmentation  technique.  The  prior 
shape  information  is  obtained  by  using  a modified  Mumford-Shah  segmentation 
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technique.  The  segmentation  is  assisted  by  a prior  shape  binary  image  and  a 
registration. 

Let  c be  a given  shape.  A binary  image  S presents  prior  shape  information 
which  is  formed  by  assigning  inside  of  c to  1 and  outside  of  c to  0.  I ia  a novel 
image  and  T is  a registration  mapping  either  a rigid  transformation  or  nonrigid 
deformation.  The  model  is  aimed  to  find  T,  ci,  and  C2  by  minimizing  the  energy 
functional: 


£(T,Cl,c2)=  [ S(x)(I(r(x))-Cl)2dx+ 

Jo. 

A /(I  - S(x))(I(  T(x))  - c2fdx  + u f |VT|2cte,  (3-12) 

Jn  J n 

where  A and  v are  positive  parameters  balancing  the  influences  from  three 
terms  in  the  model.  For  each  x € f2,  the  first  term  is  forcing  /(T(f))  to  be  close  to 
ci,  since  S(x)  = 1 for  x inside  of  c.  In  a similar  way,  the  second  term  is  compelling 
I(Y(x))  to  be  close  to  c2,  since  1 — 5(x)  = 1 for  x outside  of  c.  Smoothing  the 
vector  Y(x)  is  in  the  last  term.  After  minimizing  these  three  terms,  the  best  T, 

Ci,  and  c2  are  obtained.  In  our  numerical  experiments,  only  rigid  transformations 
are  considered,  but  the  model  can  be  generalized  to  non-rigid  deformation  as  well. 
Therefore  the  rigid  transformation  vector  T(f)  = /i Rx  + T,  where  /i  is  a scaling, 

R is  a rotation  matrix  with  respect  to  an  angle  9 , and  T is  a translation.  Since  the 
rigid  transformation  is  considered,  v = 0 in  our  experiments. 

3.2.2. 1 Euler-Lagrange  Equations  of  the  Proposed  Model 

The  evolution  equations  associated  with  the  Euler-Lagrange  equations  in 
Equation  (3-12)  are 

^ = j S(x)(I( T(x))  - ci)(VI(T(x)))Rxdx+ 

A /(I  - S(x))(/( T(x))  - c2)(VI{Y{x)))Rxdx, 

Jq 


(3-13) 
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= J S{x)(I{ T(x))  - cl){VI{T(x)))fi^xdx+ 

A j{  1 - S(x))(7( T(x))  - c2)(V7(T(x)))//^xdx,  (3-14) 


= jf  5(x)(/(T(x))  - ci)(V7(T(x)))dx+ 

A /(I  - S(x))(7(T(x))  - c2)(V7(Y(x)))dx,  (3-15) 


^ - J s(x)(I{T{x))  - cjds. 


(3-16) 


^ = - A jT(l  - 5(x))(/(T(x))  - c2)dx,  (3-17) 

where  7?  is  the  rotation  matrix  in  terms  of  the  angle  6. 

3.2.3  Numerical  Methods 

In  this  section,  numerical  methods  to  solve  the  Equations  (3-13)-(3-17)  are 
explained.  (3-12)  was  solved  by  finding  a steady  state  solution  of  (3-13)-(3-17). 

To  discretize  the  equations  (3-13)-(3-17),  we  use  a finite  difference  scheme  and  we 
applied  the  gradient  descent  method.  Let  7(T)y  denote  the  value  of  7(Y)  at  the 
pixel  Xj  = ih  and  yj  = jk,  where  h and  k is  the  pixel  size.  Denote  7(T)(x<, y2,  tn) 
by  7(T)Jf  . The  time  derivative  7(T)f  at  ( i,j,tn ) is  approximated  by  the  forward 
difference  scheme: 


7(T)t(i,  j,tn) 


7(T)n+1  _ 7(T)n 

At 


where  At  is  the  time  step.  Moreover,  (3-13)-(3-17)  were  implemented  by  the 
following  iteration  scheme  : 


(7(T)X)"  = 


7m?i+1)J  - 7(T)u- 


(/(T)y)"  = 


7(T)r0+D  - ^ 


k 
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3.2.4  Experimental  Results 

The  goal  of  this  experiment  is  to  develop  the  prior  shape  knowledge  based 
image  segmentation  technique.  The  prior  shape  information  is  obtained  by  using 
a modified  Mumford-Shah  segmentation  technique  in  section  (3.1.2).  Throughout 
experiments.  A = 1 is  used. 

The  results  of  an  application  to  the  simulated  human  brain  MRI  images  are 
from  Figure  3-10  to  Figure  3-13.  The  data  is  obtained  from  http://www.bic.rani. 
mcgill . ca/brainweb.  Simulated  MRI  volumes  for  normal  brain  with  T1  modality 
are  obtained.  Slice  thickness  is  1mm  with  9%  noise  and  20%  non-uniformity 
intensity  level.  The  goal  of  the  experiment  is  to  segment  the  light  blue  color  of 
the  simulated  brain  image.  In  Figure  3-10,  the  left  one  is  the  simulated  MRI  brain 
image  BI  and  the  right  one  is  given  binary  image  as  feature  information  S.  S is 
formed  by  assigning  the  inside  of  the  given  shape  c to  1 and  outside  0.  A novel 
image  I is  created  by  rotating  BI  by  0.1745.  In  Figure  3-11,  the  left  one  is  a novel 
simulated  MRI  brain  image  I and  the  middle  one  is  segmented  image  as  a binary 
image.  The  middle  image  in  Figure  3-11  is  created  by  applying  T to  the  prior 
shape  binary  image  S.  Finally,  the  right  one  in  Figure  3-11  is  the  boundary  of  the 
segmented  image.  After  90  iterations,  9 converges  to  -0.1763.  Initial  value  of  c\  is 
0.1  and  converges  to  0.1002  after  90  iterations.  In  a similar  way,  C2  converges  to 
0.4982  from  an  initial  value  0.5  after  90  iterations. 

Another  numerical  results  with  an  applications  to  simulated  human  brain 
image  are  shown  from  Figure  3-12  to  Figure  3-13.  The  goal  of  the  experiment  is  to 
segment  the  light  yellow  color  of  the  simulated  brain  image.  Figure  3-12,  the  left 
one  is  the  simulated  MRI  brain  image  BI  and  the  right  one  is  given  binary  image 
as  feature  information  S.  S is  formed  by  assigning  the  inside  of  the  given  shape  c 
to  1 and  outside  0.  A novel  image  I is  created  by  rotating  BI  by  -0.2618.  In  Figure 
3-13,  the  left  one  is  a novel  simulated  MRI  brain  image  I and  the  middle  one  is 
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Figure  3-10:  Left  : Simulated  MRI  image  BI,  Right  : S 


Figure  3-11:  Left  : A Novel  Image  /,  Middle  : Segmented  Image,  Right  : The 
Boundary  of  the  Segmented  Image 

segmented  image  as  a binary  image.  The  middle  image  in  Figure  3-13  is  created 
by  applying  T to  the  prior  shape  binary  image  S.  Finally,  the  right  one  in  Figure 
3-13  is  the  boundary  of  the  segmented  image.  After  155  iterations,  9 converges  to 
0.2531.  Initial  value  of  Ci  is  0.09  and  converges  to  0.0901  after  155  iterations.  In  a 
similar  way,  c2  converges  to  0.0108  from  an  initial  value  0.01  after  155  iterations. 
The  numerical  results  show  the  preliminary  effectiveness  of  the  proposed  model  to 
segment  given  image  using  a prior  shape  binary  image. 


Figure  3-12:  Left  : Simulated  MRI  image  BI,  Right  : S 
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Figure  3-13:  Left  : A Novel  Image  I,  Middle  : Segmented  Image,  Right  : The 
Boundary  of  the  Segmented  Image 

3.3  Simultaneous  Segmentation  and  Registration  using  a Modified 

Mumford-Shah  Model 

3.3.1  Introduction 

Extracting  the  boundary  and  register  the  given  images  is  one  of  the  most 
important  tasks  in  image  processing,  image  analysis,  and  computer  vision.  There 
have  been  numerous  techniques  developed  to  solve  image  segmentation  and 
registration  problems.  In  the  past,  solutions  of  these  two  problems  have  been 
studied  separately. 

Segmentation  techniques  have  been  developed  to  capture  the  object  boundary 
by  several  different  approaches;  edge-based  methods  mainly  using  active  contour 
models,  region-based  methods,  or  the  combination  of  the  two  by  using  Geodesic 
Active  Region  models.  The  most  celebrated  region  based  image  segmentation 
model  is  introduced  by  Mumford  and  Shah  in  1989  [57].  In  this  model,  an  image  is 
decomposed  into  a set  of  regions  within  the  bounded  open  set  Q and  these  regions 
are  separated  by  smooth  edges  F.  Chan  and  Vese  [9]  reformulate  a piecewise 
constant  Mumford-Shah  model  in  the  context  of  region  based  active  contour.  Chan 
and  Vese  also  proposed  a piecewise  constant  Mumford-Shah  model  in  [9,  10]  by 
using  a level  set  formulation  presented  by  Osher  and  Fedkiw  [58].  Recently,  fuzzy 
segmentation  algorithm  using  Mumford-Shah  model  is  proposed  by  Shen  [70]  and 
Thiruvenkadam  [77] . A fuzzy  segmentation  in  our  model  is  similar  to  Shen  [70] 
and  Thiruvenkadam  [77],  but  the  length  of  boundary  term  of  proposed  model  is 
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different  from  Thiruvenkadam  [77].  The  length  of  boundary  term  is  motivated  by 
Esedoglu,  et  al.  [22],  but  our  model  is  variational  which  differs  from  Esedoglu,  et 

al-  [22). 

Image  registration  is  the  process  of  overlaying  two  or  more  images  of  the  same 
scene  taken  at  different  times,  and/or  by  different  sensors.  Area  based  methods  or 
feature  based  methods  have  been  developed  to  match  given  images.  For  the  details 
of  the  image  registration  techniques,  refer  to  [92]. 

A joint  segmentation  and  registration  methods  have  been  developed  by  Chen, 
Yezzi.  Young,  et  al.  [2,  13,  15,  88,  89,  90].  An  explicit  combination  of  registration 
with  segmentation  has  been  developed  in  a variational  framework  through  active 
contours  by  Yezzi.  Zollei,  and  Kapur  [88,  89].  Moelich  and  Chan  [56]  extend 
the  algorithms  of  Yezzi,  Zollei,  and  Kapur  [88]  to  joint  segmentation  and  of  an 
object  in  two  images,  and  the  prior  knowledge  is  implicitly  incorporated  through 
the  simultaneous  evolution  of  the  registration  maps  and  the  active  contour  in 
the  variational  framework.  The  morphing  active  contours  algorithm  is  combined 
with  the  joint  segmentation  and  registration  with  an  application  to  CT  scans  for 
radiotherapy  treatment  planning  presented  by  Levy,  et  al.  [90]. 

The  purpose  of  the  model  is  segment  and  register  given  images  simultaneously 
using  a modified  Mumford-Shah  technique  and  region  intensity  values.  The  model 
performs  a simultaneous  segmentation  and  registration  in  a similar  way  from  Yezzi, 
et  al.  [88,  89].  But  the  model  differs  from  Yezzi,  et  al.  [88,  89]  by  using  a region 
intensity  and  a modified  Mumford-Shah  segmentation  technique.  This  section  is 
organized  as  follows:  In  section  3.3.2,  the  model  is  proposed  for  a simultaneous 
segmentation  and  registration  and  the  Euler- Lagrange  equations  of  the  proposed 
model  is  described.  In  section  3.3.3,  numerical  methods  pertaining  to  the  model 
are  explained.  Experimental  results  of  the  model  which  were  applied  to  synthetic 
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and  human  brain  MRI  images  are  also  showed  in  this  section.  In  section  3.3.4,  the 
conclusion  follows  and  future  work  is  stated. 

3.3.2  Description  of  the  Proposed  Model 

A new  variational  region  based  model  for  a simultaneous  image  segmentation 
and  registration  using  a modified  Mumford-Shah  technique  is  suggested.  The 
purpose  of  the  model  is  segment  and  register  given  images  simultaneously  using  a 
modified  Mumford-Shah  technique  and  region  intensity  values.  The  segmentation 
is  obtained  by  minimizing  a modified  Mumford-Shah  model.  A global  rigid 
registration  is  assisted  by  the  segmentation  information  and  region  intensity  value. 
A segmentation  and  registration  process  is  obtained  simultaneously  in  this  model. 
The  model  is  aimed  to  find  0,  ci,  C2,  T,  dx,  and  d2  by  minimizing  the  energy 
functional: 

£’(0,c1,c2,T,d1,d2)  = Ai  [ (sin2Q(x)(Ix  ~ c\)2  + cos2Q(x)(Ix  - c2)2)dx+ 

Jn 

f («m’( 29(x))| vei2  + We(x)co»<e(j:))|a  + xf , vr|2(te 

Jn  e Jn 

A4  f {sin2Q(x)(I2( T(x))  - dx)2  + cos2Q(x)(I2(Y(x ))  - d2)2)dx  (3-18) 

Jn 

where  Ix  and  /2  are  novel  images,  Q is  domain,  and  A,  > 0 (i  = 1,2, 3, 4)  are 
parameters  balancing  the  influences  from  the  four  terms  in  the  model. 

The  first  term  forces  sin2(0(i)),  towards  0 if  h(x)  is  different  from  a and 
towards  1 if  I\{x)  is  close  to  Ci,  for  every  x G fl.  For  .r  G Q,  the  average  intensity 
of  I\{x)  is  close  to  c1(  where  sin2(Q(x))  > cos2(Q(x)).  In  a similar  way,  cos2(0(i)), 
towards  0 if  I\ (x)  is  different  from  c2  and  towards  1 if  I\{x)  is  close  to  c2,  for  every 
x E Q.  For  x G fl,  the  average  intensity  of  Ix(x)  is  close  to  c2,  where  sin2(Q(x))  < 
cos~(Q(x)).  The  second  term  is  a length  of  boundary  term.  Smoothing  the  vector 
T(x)  is  in  the  third  term.  For  each  x G fl,  the  fourth  term  is  forcing  sin2(Q(x)), 
towards  0 if  72(T(x))  is  different  from  dx  and  towards  1 if  72(T(x))  is  close  to  dx. 
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For  x E El,  the  average  intensity  of  I2(T(x))  is  close  to  dx,  where  sin2(Q(x))  > 
cos2(Q(x)).  In  a similar  way,  cos2(6(x)),  towards  0 if  I2( T(x))  is  different  from 
d2  and  towards  1 if  /2(T(x))  is  close  to  d2.  For  x E El,  the  average  intensity  of 
72(Y(a:))  is  close  to  d2 , where  sm2(0(x))  < cos2(0(x)).  After  minimizing  these 
three  terms,  the  best  0,  T,  c\,  c2,  d\ , and  d2  are  obtained.  For  every  x E El,  denote 
<t>{x)  = sm2(0(i))  — cos2(Q(x )).  If  <p  > 0,  then  assign  1 to  inside  of  the  given  shape 

c.  If  0 < 0.  then  assign  0 to  inside  of  the  given  shape  c.  The  segmented  binary 
image  of  I\  is  formed.  In  a similar  way,  H (</>(¥))  segments  a novel  image  I2  by 
forming  a binary  image  which  has  a 1 inside  of  the  given  shape  d and  0 outside  of 

d. 

In  our  numerical  experiments,  only  rigid  transformation  is  considered,  but 
the  model  can  be  generalized  to  non-rigid  deformation  as  well.  Therefore  the  rigid 
transformation  vector  T(x)  = /iRx  + T , where  //  is  a scaling,  R is  a rotation  matrix 
with  respect  to  an  angle  9 , and  T is  a translation.  Since  the  rigid  transformation  is 
considered,  A3  = 0 in  our  experiments. 

3.3.2. 1 Euler-Lagrange  Equations  of  the  Proposed  Model 

The  evolution  equations  associated  with  the  Euler-Lagrange  equations  in 
Equation  (3-18)  are 


— = Ai(sm(20(x))(/i  - ci)2  - sin(2Q(x))(Ii  - c2)2)dx 

-{2esm(40(£))|V0|2  - 2esm2(20(j))A0 
+ A24sm30(x)cos30(x)cos(20(x)) 


+A4(sm(20(x))(/2  - di)2  - sin(2Q(x))(I2  - d2)2)dx,  in  El 

<90 

——  = 0,  on  dEl 


(3-19) 


(3-20) 
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dc  F 

= — Ai  J cos20(x)(/1( T(x))  - c2)dx. 

— A4  J (sin(2G(x))(I2(Y(x)))  - di)2  - sin(2Q(x)) 

(I2(Y(x)))  - d2)2)Rxdx  + A4  f sin2Q(x) 

J n 

(/2(T(x))  - di)(VJa(T(*)))J2x+ 
cos20(x)(/2(T(x))  - d2)(V/2(T(x)))/&dx, 


cd 


A4  f (sm(20(x))(/2(T(x)))  — di)2  - sm(20(x)) 
J n 


(/2(T(x)))  - d2)2)fi^xdx  + A4  [ sin2Q(x) 
dv  J q 

(/2(T(x))  - d1)(V/2(T(x)))/2^x+ 

H R 

cos2G(x)(I2( T(x))  - d2)(V/2(T(x)))^— xdx. 


nm  /» 

= A4  J (sin(20(x))(/2(Y(x)))  - d])2  - sm(29(x)) 

(/2(T(x)))  - d2)2)dx  + A4  f sin2 Q(x) 

Jfi 

(I2(Y(x))  - d1)(VI2(Y(x)))+ 

cos2G(x)(/2(Y(x))  - d2)(V/2(T(x)))dx, 


dd  F 

= -A4  J sin2Q(x)(I2(Y(x))  - dx)dx, 
dd  f 

= — A4  J cos2Q(x){I2( T(x))  - d2)dx, 
where  /?  is  the  rotation  matrix  in  terms  of  the  angle  d. 


(3-21) 

(3-22) 

(3-23) 

(3-24) 

(3-25) 


(3-26) 
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3.3.3  Numerical  Methods 

In  this  section,  numerical  methods  to  solve  the  Equations  (3-19)-(3-26)  are 
explained.  (3-18)  was  solved  by  finding  a steady  state  solution  of  (3-19)-(3-26). 

To  discretize  the  equations  (3-19)-(3-26),  we  use  a finite  difference  scheme  and  we 
applied  the  gradient  descent  method.  Let  0SJ  denote  the  value  of  0 at  the  pixel 
= ih  and  y,  = jk,  where  h and  k is  the  pixel  size.  Denote  9(xi,yj,tn)  by  0^.  The 
time  derivative  0,  at  ( i,j,tn ) is  approximated  by  the  forward  difference  scheme: 

0 ft  it)-  ^lJ 

Ji  t n)  a , t 


where  At  is  the  time  step.  /(Y)jj  denote  the  value  of  /( T)  at  the  pixel  = ih  and 
Uj  — jk , where  h and  k is  the  pixel  size.  Denote  /2(Y)(:ri,  y,,  t„)  by  72(T)?.  The 
time  derivative  72(Y)f  at  ( i,j,tn ) is  approximated  by  the  forward  difference  scheme: 


/2(T  )t(i,j,tn)  = 


frcny1  - 72(T)5 
At 


where  At  is  the  time  step.  Moreover,  (3-19)-(3-26)  were  implemented  by  the 
following  iteration  scheme  : 


(e.)&  = 


fln  an 


(e„)S  = 


(~yt  _ on 
Wt(j+1)  uu 


,J  k 

n \n  ©IUlM-20S+0?i.  1M 

pin  opin  _l  ftn 

to.  \n  _ ” (j+D  + ”0-1) 

x^wJij  ~ fc2 

(«T),)'‘  = /z(T)t'A-1)j  - 

(/,(TW8  - 


Zb 
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3.3.4  Experimental  Results 

The  purpose  of  this  experiment  is  segment  and  register  given  images  simul- 
taneously using  a modified  Mumford-Shah  technique  and  region  intensity  values. 
The  segmentation  is  obtained  by  minimizing  a modified  Mumford-Shah  model. 
Throughout  the  experiments,  A:  = 3,  A2  = 1/9,  and  A4  = 1 are  used.  The  results  of 
an  application  to  the  simulated  human  brain  MRI  images  are  from  Figure  3-14  to 
Figure  3-17.  The  data  is  obtained 

from  http://www.bic.mm.mcgill.ca/brainweb.  Simulated  MRI  volumes  for 
normal  brain  with  T1  modality  are  obtained.  Slice  thickness  is  1mm  with  9% 
noise  and  20%  non-uniformity  intensity  level.  For  Figure  3-14  and  Figure  3-15,  a 
given  image  I2  is  created  by  rotating  an  image  /x,  by  amount  of  0.1745.  In  Figure 
3-14.  from  left  to  right,  the  first  one  is  a novel  image  7i,  the  second  one  is  a novel 
image  I2,  and  the  third  one  is  a segmented  image  of  I\.  For  every  x 6 f2,  denote 
< t>{x ) = sin2(Q(x))  — cos2(Q(x)).  If  4>  > 0,  then  assign  1 to  inside  of  the  given 
shape  c.  If  </>  < 0,  then  assign  0 to  inside  of  the  given  shape  c.  The  segmented 
binary  image  of  Ii  is  formed  by  H(cf)(x)).  In  Figure  3-15,  from  left  to  right,  the  first 
one  is  the  boundary  of  the  segmented  image  of  . The  second  one  is  segmented 
binary  image  of  /2  formed  by  //(</>( T)).  Finally,  the  third  one  is  the  boundary  of 
the  segmented  image  of  /2.  After  20  iterations,  6 converges  to  -0.1776.  Initial  value 
of  ci  is  0.01  which  converges  to  0.1173  after  20  iterations  and  c2  converges  to  0.5907 
from  an  initial  value  0.04  after  20  iterations.  In  a similar  way,  initial  value  of  d\  is 
0.02  which  converges  to  0.1404  after  20  iterations  and  d2  converges  to  0.5393  from 
an  initial  value  0.05  after  20  iterations. 

Finally,  a given  image  I2  is  created  by  rotating  an  image  /x , by  amount  of 
-0.2618  for  Figure  3-16  and  Figure  3-17.  In  Figure  3-16,  from  left  to  right,  the  first 
one  is  a novel  image  Ix,  the  second  one  is  a novel  image  /2,  and  the  third  one  is  a 
segmented  image  of  I\.  For  every  x e Cl,  denote  <f>(x)  = sin2(Q(x ))  - cos2(0(x)).  If 
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Figure  3-14:  Left  : A Novel  Image  A,  Middle  : A Novel  Image  I2,  Right  : Seg- 
mented Image  of  I\ 


Figure  3-15:  Left  : The  Boundary  of  the  Segmented  Image  A,  Middle  : Segmented 
Image  of  /2,  Right  : The  Boundary  of  the  Segmented  Image  of  I2 

4>  > 0,  then  assign  1 to  inside  of  the  given  shape  c.  If  cj)  < 0,  then  assign  0 to  inside 
of  the  given  shape  c.  The  segmented  binary  image  of  I\  is  formed  by  H((f>(x)).  In 
Figure  3-17,  from  left  to  right,  the  first  one  is  the  boundary  of  the  segmented  image 
of  I\.  The  second  one  is  segmented  binary  image  of  I 2 formed  by  H((f){ T)).  Finally, 
the  third  one  is  the  boundary  of  the  segmented  image  of  I2.  After  20  iterations,  9 
converges  to  0.2568.  Initial  value  of  Ci  is  0.45  which  converges  to  0.4496  after  20 
iterations  and  c2  converges  to  0.0293  from  an  initial  value  0.2  after  20  iterations. 

In  a similar  way,  initial  value  of  d\  is  0.55  which  converges  to  0.4263  after  20 
iterations  and  d2  converges  to  0.0548  from  an  initial  value  0.4  after  20  iterations. 
The  experimental  results  show  the  effectiveness  of  the  model  in  detecting  the 
boundaries  of  the  given  objects  and  registering  novel  images  simultaneously  and 
robustness  to  a noise. 
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Figure  3-16:  Left  : A Novel  Image  I\,  Middle  : A Novel  Image  1 2,  Right  : Seg- 
mented Image  of  I\ 


Figure  3-17:  Left  : The  Boundary  of  the  Segmented  Image  of  /1;  Middle  : Seg- 
mented Image  of  /2,  Right  : The  Boundary  of  the  Segmented  Image  of  J2 


CHAPTER  4 

CONCLUSIONS  AND  FUTURE  DIRECTION 
In  Chapter  2,  a new  variational  method  for  a simultaneous  segmentation  and 
non-rigid  registration  that  incorporates  the  shape  and  intensity  priors  is  proposed. 
Since  the  heart  cardiac  images  can  be  varied  images  to  images,  a local  non-rigid 
deformation  term  is  added  to  capture  the  boundary  effectively.  Numerical  results 
showed  the  effectiveness  of  the  suggested  model. 

The  prior  shape  C*  and  prior  image  I*  which  were  created  by  using  the 
method  from  [12]  are  in  Figure  2-1.  The  average  contour  is  used  as  a prior  shape 
and  an  average  image  for  a prior  intensity.  From  Figure  2-2  to  Figure  2-3,  the 
effectiveness  of  each  term  in  the  model.  From  Table  2-1,  it  is  easily  seen  that  the 
major  contribution  terms  in  the  model  are  the  first  and  third  terms.  From  Figure 
2-4  to  Figure  2-7,  the  model  performance  to  the  given  images  which  is  compared  to 
expert  borders  is  shown.  The  global  rigid  transformation  model  was  calculated  in 

u 


Equation  (2.9)  without 


term.  The  distance  measurements  of  the  results  of 


the  model  from  the  expert’s  borders  by  using  the  distance  function  are  followed  in 
Table  2-2  and  Table  2-3.  The  results  show  the  effectiveness  of  the  model  to  capture 
the  boundary  of  the  given  image,  which  is  better  than  the  model  with  global  rigid 
transformation  term  only.  These  results  were  obtained  through  Matlab  and  the 
computational  time  in  each  iteration  was  about  1 second  using  a Pentium  4 CPU 
running  at  2.4  GHZ  with  512  MB  of  RAM.  Windows  XP  Home  Edition  was  used 
as  the  operating  system. 

Due  to  the  inhomogeneity,  dropout,  and  loss  of  some  information  near  the 
edges  in  heart  images,  parameters  are  adjusted  in  each  image  to  get  desired  results. 
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Some  of  the  images  with  a diminishing  problem  which  depends  on  the  parameters 
are  observed  in  the  experiments.  Since  the  shape  prior  is  used  as  an  initial  contour, 
the  level  set  form  of  the  model  is  fixed  which  does  not  require  the  re-initialization 
process  during  the  numerical  calculation.  The  intensity  matching  term  in  the 
energy  function  of  the  model  provides  the  difference  of  the  intensity  variation 
across  the  prior  shape  and  the  shape  of  interest.  Minimizing  the  difference  of  the 
gradients  of  the  intensity  between  smoothed  prior  image  and  the  smoothed  given 
image  is  more  accurate  than  minimizing  the  gradient  term  in  the  given  image,  due 
to  the  loss  of  the  information  of  the  gradient  near  the  edge  of  the  desired  object. 
Moreover,  this  model  can  be  extended  to  3D  cases  and  any  other  types  of  images. 
Generalization  of  the  model  to  global  non-rigid  deformation  is  needed  in  future 
work. 

Several  statistical  algorithms  for  generating  shape  models  are  presented.  The 
aim  of  the  models  is  to  develop  an  approach  for  the  mean  shape  and  clustering 
to  detect  differences  in  the  shape  of  anatomical  structures.  One  of  the  suggested 
models  minimizes  a shape  related  energy  to  create  the  mean  shape  with  the 
orthogonal  Procrustes  method.  The  model  does  not  depend  upon  the  choice  of 
the  initial  contour  and  generates  the  mean  shape  and  alignments  to  the  mean 
shape  from  the  given  contours  simultaneously.  The  other  models  generate  the  mean 
shape  and  clustering  by  using  Self-Organizing  map.  The  closeness  is  measured 
by  the  orthogonal  Procrustes  method  or  area  difference  distance.  This  technique 
has  the  advantage  of  generating  the  mean  shape  of  each  cluster  and  clustering 
given  contours  simultaneously.  In  addition,  Self-Organizing  map  combined  with  a 
Principal  Component  Analysis  (PC A)  model  is  also  introduced.  Numerical  results 
with  an  application  to  human-heart  normal  two-chamber  end-diastolic  epicardial 
borders  are  shown  the  effectiveness  of  the  suggested  models.  Figure  2-8  represents 
the  shape  related  energy  function  model  result  applied  eighty-five  human-heart 
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end-diastolic  epicardial  borders.  For  Figure  2-8,  a dotted  line  in  the  left  one  is 
the  initial  contour  of  the  given  contours.  The  dotted  line  in  the  right  one  is  mean 
shape  of  the  given  contours.  Three  distance  measurements,  namely  mean  distance, 
standard  deviation,  and  Tanimoto  measurement  of  the  result  of  the  models  from 
Figure  2-8  is  in  Table  2-6.  From  Figure  2-9  to  Figure  2-13  are  the  results  which 
were  obtained  when  the  Self-Organizing  map  was  applied  to  the  eighty-five  human- 
heart  two-chamber  normal  end-diastolic  epicardiac  borders  to  obtain  clustering 
and  the  mean  shape  of  each  cluster.  Figure  2-9  is  the  numerical  results  using 
Self-Organizing  map  with  the  area  difference  distance  method.  In  a similar  way, 
Figure  2-10  is  the  numerical  results  using  Self-Organizing  map  with  orthogonal 
Procrustes  method.  It  is  easily  seen  that  the  orthogonal  method  is  more  efficient 
than  area  difference  method  from  Table  2-7.  To  get  better  mean  shape,  Figure  2-11 
shows  the  numerical  results  using  the  method  from  section  2.8.2. 3.  Pre-alignments 
and  applying  principal  component  analysis  give  better  results  which  can  be  easily 
seen  in  Table  2-7.  From  Figure  2-12  and  Figure  2-13,  numerical  results  with  an 
application  to  the  eighty-five  human-heart  two-chamber  normal  end-diastolic 
epicardiac  borders  to  cluster  into  three  groups.  Figure  2-12  shows  the  results  using 
the  method  from  section  2. 8. 2. 2 with  orthogonal  Procrustes  method.  First  group, 
there  are  30  contours,  28  contours  are  in  the  second  group.  Finally  27  contours  are 
in  the  third  group.  Numerical  results  using  the  method  from  section  2. 8. 2. 3 is  in 
Figure  2-13.  Distance  measurements  comparisons  from  Table  2-8  and  Table  2-9 
show  the  effectiveness  of  the  suggested  models. 

Chapter  3 contains  three  interrelated  sections;  Image  segmentation  using  a 
modified  Mumford-Shah  technique,  Region  based  image  segmentation  with  prior 
shape,  and  a modified  Mumford-Shah  model  based  simultaneous  segmentation  and 
registration. 
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A region  based  variational  partial  differential  equation  (PDE)  model  for  an 
image  segmentation  is  presented  with  an  application  to  synthetic  and  human  brain 
MRI  images.  The  purpose  of  the  model  is  to  segment  borders  of  the  given  images. 
The  segmentation  is  obtained  by  using  a modified  Mumford-Shah  segmentation 
technique.  The  numerical  results  with  an  application  to  the  synthetic  image  and 
simulated  normal  human  brain  images  are  shown  from  Figure  3-1  to  Figure  3-9. 
Preliminary  numerical  results  show  the  effectiveness  of  the  model  to  capture  the 
boundary  of  the  given  fuzzy  images. 

A new  variational  region  based  image  segmentation  model  based  on  the  region 
intensity  and  prior  shape  information  is  suggested.  The  prior  shape  information  is 
extracted  by  using  a modified  Mumford-Shah  segmentation.  The  shape  knowledge 
information  and  registration  are  supporting  the  image  segmentation  process.  The 
model  results  with  an  application  to  simulated  human  brain  images  are  shown  from 
Figure  3-10  to  Figure  3-13.  These  numerical  results  show  the  effectiveness  of  the 
suggested  model  to  register  given  wrongly  mapped  images  correctly.  Even  though 
numerical  experiments  are  done  only  for  the  global  rigid  transformation,  it  can  be 
generalized  to  global  non-rigid  deformation.  Generalization  of  the  model  to  global 
non-rigid  deformation  is  needed  for  future  work. 

A new  variational  region  based  model  for  a simultaneous  image  segmentation 
and  registration  using  a modified  Mumford-Shah  technique  is  presented.  The 
purpose  of  the  model  is  segment  and  register  given  images  simultaneously  using  a 
modified  Mumford-Shah  technique  and  region  intensity  values.  The  model  results 
with  an  application  to  simulated  human  brain  images  are  shown  from  Figure 
3-14  to  Figure  3-17.  The  numerical  results  show  the  effectiveness  of  the  model 
in  detecting  the  boundaries  of  the  given  objects  and  registering  novel  images 
simultaneously.  This  model  can  be  extended  to  3D  cases  and  any  other  types  of 


images. 
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The  discovery  of  X-rays  by  German  physicist  Rontgen  in  1895  initiated 
the  history  of  medical  imaging.  Since  then,  many  medical  devices  have  been 
developed  to  diagnose  and  treat  different  diseases  or  to  aid  doctors  during  surgery. 
Image  qualities  which  are  produced  by  medical  devices  are  not  good  enough 
to  be  applied  to  clinical  research,  therefore  the  development  of  mathematical 
algorithms  has  become  necessary.  The  purpose  of  this  dissertation  is  to  discover 
the  mathematical  algorithms  needed  to  improve  these  imaging  problems,  especially 
with  an  application  to  human  heart  and  brain  images.  Suggested  models  in  this 
dissertation  with  applications  to  medical  imaging  data  show  some  promising 
results.  In  addition,  these  algorithms  can  be  applied  not  only  to  medical  imaging 
but  also  to  general  imaging. 


